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EVALUATION  OF  DENSITIES  AND  DISTRIBUTIONS  VIA  HERMITE  AND  GENERALIZED 
LAGUERRE  SERIES  EMPLOYING  HIGH-ORDER  EXPANSION  COEFFICIENTS 
DETERMINED  RECURSIVELY  VIA  MOMENTS  OR  CUMULANTS 

INTRODUCTION 


'In  the  theoretical  analysis  of  performance  of  some  systems  with 
nonlinearities  and/or  memory,  it  often  happens  that  the  only  statistics  about 
the  decision  (or  output)  random  variable  of  interest  that  can  be  easily  found 
are  the  moments,  or  in  other  cases,  the  cumulants.  Explicit  relations  for  the 
low-order  expansion  coefficients  in  Edgeworth  or  Gram-Charl ier  series  are 
available  in  terms  of  the  available  moments  or  cumulants  [1,  pp.  172  and  191], 
[2,  pp.  223  and  226],  [3,  pp.  157  and  159].  However,  for  higher-order  moments 
and  cumulants,  these  explicit  nonrecursive  relations  are  very  tedious  to 
derive,  become  extremely  lengthy,  and  are  not  practical  to  use. 

We  will  address  the  problem  of  obtaining  accurate  high-order  series 
expansion  approximations  of  the  probability  density  function  and  cumulative 
distribution  function  of  a  random  variable  of  interest,  in  terms  of  the 
available  moments  or  cumulants  of  that  random  variable.  The  necessity  of 
being  able  to  approximate  probability  density  functions  and  cumulative 
distribution  functions  from  knowledge  of  either  the  moments  or  the  cumulants, 
is  that  some  physical  problems  have  these  particular  statistics  as  natural  and 
convenient  starting  points.  For  example,  if  a  physical  processor  sums 
together  a  number  of  independent  Rician  random  variates,  the  characteristic 
function  and  cumulants  of  the  individual  random  variables  or  their  sum  are  not 
available  in  any  useful  analytic  form;  however,  the  high-order  moments  of  an 
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individual  Rician  variate  can  be  easily  and  accurately  evaluated  by 
recurrence,  and  thereby  the  moments  of  the  sum  can  be  obtained.  Conversely, 
for  shot  noise  with  random  amplitude  and  duration  modulation,  the  probability 
density  function  is  not  readily  available,  whereas  the  characteristic  function 
is,  and  the  cumulants  are  simple  to  evaluate  [4,  appendix  C]. 

The  particular  series  expansions  we  employ  are  based  on  the  two  special 
classes  of  weighting  functions 

w(u)  = — ^-rr—  exp/-  for  all  u  Hermite  (1) 

(2ir  {  2b  ) 

and 

w(u)  =  — — \illl  for  u  >  o  generalized  Laguerre.  (2) 

Ba  1  Fun) 

The  orthonorrnal  polynomials  associated  with  these  weightings  are  directly 
related  to  the  Hermite  and  generalized  Laguerre  polynomials,  respectively 
[5,  22.2.15  and  22.2.13].  The  weightings  each  have  two  free  parameters,  a  and 
B,  which  can  be  manipulated  to  advantage  in  obtaining  finite  (high-order) 
series  expansions  which  well  approximate  a  given  (unknown)  probability  density 
function  and  cumulative  distribution  function. 

The  question  of  when  a  set  of  moments  uniquely  determines  the  probability 
density  function  is  a  difficult  one;  see,  for  example,  [3,  pp.  109-112  and 
179].  Also,  the  convergence  of  the  series  is  very  involved  [2,  pp.  223  and 
258],  [3,  pp.  161-163].  But,  even  if  the  series  is  divergent,  use  of  a 
limited  number  of  expansion  coefficients  often  gives  a  satisfactory 
approximation  to  the  desired  probability  density  function  [3,  p.  167].  We 
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presume  here  that  the  moments  do  uniquely  determine  the  probability  density 
function  and  are  buoyed  in  that  respect  by  the  comment  [3,  p.  87]  that  most 
distributions  in  statistical  practice  do  possess  this  property. 

The  main  idea  in  the  series  expansion  approach  here  is  not  necessarily  to 
get  as  many  terms  as  possible,  but  rather  to  get  as  rapid  convergence  as 
possible  of  the  series.  If  a  particular  choice  of  weighting  parameters  a  and 
6  results  in  sufficiently  small  expansion  coefficients,  say,  at  order  10,  this 
is  better  than  another  choice  of  a  and  s  where  20  or  30  terms  are  required  for 
the  same  size  coefficients.  In  fact,  if  a  and  e  could  be  chosen  such  that  the 
series  terminated  (zero  coefficients)  after  a  few  terms,  that  would  be  ideal; 
however,  this  is  not  the  case,  and  in  fact,  the  choice  of  a  and  @  requires 
some  trial-and-error  to  achieve  rapidly  decreasing  coefficients. 

The  expansion  coefficients  of  a  given  probability  density  function,  in  an 
orthonormal  set  of  Hermite  or  generalized  Laguerre  polynomials,  are  denoted  by 
l^n^O’  w^ere  N  is  the  number  of  available  or  known  moments  or  cumulants. 

Very  often,  the  choice  of  a  and  3  in  (1)  or  (2)  has  been  made  such  that 
b,  =  0  and  b2  =  0,  for  purposes  of  analytic  simplicity  and  for  hopeful 
early  termination  of  the  series;  see  for  example  [1,  pp.  171  and  191], 

[2,  p.  223],  [3,  p.  159].  However,  it  will  be  demonstrated  that  this  is 
generally  not  the  best  choice,  and  that  more  rapidly  decaying  coefficients  can 
be  achieved  by  other  (mismatched)  values  of  a  and  3,  which  must  be  searched 
for  numerically;  this  possibility  is  also  mentioned  in  [3,  p.  164].  In  fact, 
an  example  will  be  given  which  illustrates  that  the  choice  of  parameters  a  and 
B  to  make  expansion  coefficients  b^^  and  b2  zero,  can  in  fact,  lead  to  a 
divergent  Hermite  series. 
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Depending  on  the  available  information  about  the  probability  density 
function,  i.e.,  moments  or  cumulants,  a  variety  of  methods  will  be  given  for 
determining  the  expansion  coefficients  [b^.  In  particular,  for  both  the 
Hermite  and  generalized  Laguerre  series,  we  can  get  the  coefficients  by  three 
different  procedures: 

(a)  recursively  via  cumulants, 

(b)  directly  via  moments, 

(c)  recursively  via  moments. 

The  reason  for  having  these  alternatives  is  that  the  calculation  of  expansion 
coefficients  {b^  for  high-order  n  invariably  runs  into  large  round-off 
error.  In  order  to  reduce  this  round-off  error,  the  amount  of 
number-crunching  on  the  computer  should  be  minimized,  and  any  spurious 
transformations  between  moments  and  cumulants  should  be  avoided  if  possible. 
Thus  it  is  desireable  to  have  techniques  which  can  accomplish  the  desired  goal 
of  evaluating  expansion  coefficients  jb  as  directly  as  possible  from  the 
available  information.  The  use  of  different  alternatives  also  enables 
comparisons  of  the  computed  expansion  coefficients  and  thereby  furnishes 
quantitative  assessment  of  the  amount  of  round-off  error.  Recursive 
inter-relationships  between  moments,  central  moments,  and  cumulants  are  given 
in  [6],  including  cases  of  two  dependent  random  variables. 
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FUNDAMENTAL  EQUATIONS 


DEFINITION  OF  STATISTICS 

Suppose  a  function  p  has  known  moments* 

un  =  J  du  un  p(u)  for  0  £  n  <_  N  .  (3) 

The  function  p  need  not  have  unit  area,  i.e.,  Uq  ^  1  is  allowed,  and  p  can 

become  negative  at  some  arguments  u.  Nevertheless,  for  convenience,  and  since 

most  of  our  applications  are  to  random  variables,  we  shall  refer  to  p  as  a 

probability  density  function,  and  to  its  running  integral 

a 

P(u)  =  j'  dt  p(t)  (4) 

—  O0 

as  a  cumulative  distribution  function.  We  shall  presume  that  Uq  >  0  in  all 
cases. 


The  characteristic  function  corresponding  to  probability  density  function 
p  is  the  Fourier  transform 

(5) 


(6) 


f (if )  =  du  exp(ifu)  p(u)  . 


When  f  is  expanded  in  a  power  series,  the  result  is 

<*> 


f(i?)  = 

n=0 


*  Integrals  without  limits  are  over  the  range  of  nonzero  integrand. 
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in  terms  of  the  moments  in  (3).  Alternatively,  if  Jin  f  is  expanded  in  a  power 
series, 

°o 

>  f Cl?)  -  2  %,  ■  (?) 

n=0 

where  the  quantities  are  the  cumulants  of  p  or  f.  Observe  that 
generally,  to  the  lowest  three  orders. 


XQ  =>  f(0)  =  Jin  u0  *  o 


(8) 


The  available  information  on  probability  density  function  p  will  be 
either 


moments 


or  cumul 


ants  tX^ 


N 

0  • 


(9) 


Whichever  is  available,  we  wish  to  get  high-order  accurate  approximations  to  p 
and  cumulative  distribution  function  P  in  (4);  that  is,  values  of  N  in  the 
order  of  10  to  100  are  of  interest. 
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WEIGHTING  FUNCTION  PROPERTIES 


We  select  a  nonnegative  weighting  function  w  such  that 
w(u)  >  0  at  least  where  p(u)  ^  0  . 


(10) 


We  also  disallow  any  impulses  in  w.  The  moments  of  weighting  w  are  defined 
analogously  to  (3)  as 


(11) 


for  n  >  0  ; 


it  is  presumed  that  these  quantities  can  be  evaluated  for  as  large  n  as 
required. 

Suppose  weighting  w  has  r  free  parameters  (plus  a  scaling  parameter).  It 
might  then  seem  beneficial  to  choose  them  such  that  the  moments  of  w  and  p  are 
approximately  equal. 


vn  £  un  for  1  <_  n  <  r  (plus  vQ  2  uQ)  , 


(12) 


for  then  the  abscissa  scales  of  w  and  p  would  tend  to  match.  However,  (12) 
will  turn  out  to  be  not  so  desireable,  and  the  choice  of  the  r  weighting 
parameter  values  should  be  based  on  another  criterion.  The  ordinate  scale  of 
w  is  actually  immaterial,  since  the  expansion  coefficients  {b  1  will  absorb 
this  scaling;  so  henceforth  we  presume  that  vq  =  1  with  no  loss  of 
general ity. 
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APPROXIMATION  PROCEDURE 


Let  be  any  n-th  order  polynomial,  and  approximate  probability 
density  function  p  by  function 

N 

PN(u)  =  w(u)  bn  Qn(u)  where  w(u)  >  0  , 


where  [b^Q  are  the  expansion  coefficients.  Define  weighted  squared 
error 


=  J*  du  y(u)  [p(u)-pN(u)]2  = 


J  du  Y  ( u )  [p  ( u )  — w  ( u ) 


b 

n 


(u)]2 


(13) 


(14) 


where  error-weighting  y  is  nonnegative.  If  we  minimize  Ejs,  by  choice  of 
expansion  coefficients  [b^Q,  there  follows  the  set  of  linear  equations 

b4du 

In  order  to  use  only  the  available  information  in  (9)  about  p,  the  right-hand 
side  of  (15)  must  simplify  according  to  the  selection 

y(u)  =  — ^--y  where  w(u)  >  0  (and  arbitrary  elsewhere).  (16) 

Furthermore,  since  constant  K  merely  scales  error  E^,  and  appears  on  both 
sides  of  (15),  we  can  set  K  =  1  without  loss  of  generality.  Then  (14)  becomes 


y(u)  w^(u)  Qk(u)  Qn(u)  = 


du  y(u)  w ( u )  p ( u )  Qk(u) 


(15) 

for  U  <  k  <  N. 
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‘N 


=  Idu  w(u)fe$ 


n=0 


b  Q  ( 
n  xrr 
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where  w(u)  >  0  ,  (17) 


and  (15)  reduces  to 
N 
n=i 


bn Jdu  w(u)  Qk(u)  Qn(u)  =  Jdu  p(u)  Qk(u)  for  0  £  k  £  N.  (18) 


In  general,  this  is  N+l  simultaneous  linear  equations  in  the  N+l  unknowns 
r  1  N  k 

{bnj0.  The  choice  Qk(u)  =  u  would  lead  to  an  apparently  simple  set 
of  equations,  when  (11)  and  (3)  are  used.  However,  a  few  numerical  examples 
quickly  reveals  that  they  are  very  ill-conditioned,  due  to  the  character  of 
the  nondiagonal  matrix  with  elements 

J'du  w(u)  Qk(u)  Qn(u)  for  0  £  k,n  £  N  (19) 


that  appears  on  the  left-hand  side  of  (18).  In  order  to  avoid  the  significant 
round-off  error  associated  with  solving  such  a  system  for  large  N,  we  choose 
^n\|j  to  a  set  of  orthonormal  polynomials  with  respect  to  weighting 
w;  i.e.,  (19)  is  1  for  k  =  n,  and  0  otherwise.  Also  recall  that 
v0  =  jdu  w(u)  =  1  without  loss  of  generality. 


Equation  (18)  then  reduces  to  an  explicit  relation  for  the  expansion 
coefficients : 


yu> 


for  0  <  k  <  N  , 


(20) 


and  (17)  for  the  error  becomes  merely 


du 


P2(u) 

w(u) 


b 


2 

n 


(21) 
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It  will  be  presumed  that  the  integral  in  (21)  is  finite;  otherwise,  the 
error  would  be  infinite,  which  is  a  meaningless  problem.  This  will  put  some 
restrictions  on  the  parameter  choices  of  weighting  w,  since  this  error 
integral  depends  on  these  parameters  as  well  as  on  the  given  probability 
density  function  p.  The  sum  of  squares  in  (21)  must  then  be  bounded,  and  in 
fact  affords  a  measure  of  the  adequacy  of  approximation  (13),  by  saturating 
(at  an  apriori  unknown  value)  for  large  N. 

As  N  increases,  the  values  of  the  lower-order  expansion  coefficients 
|b^  in  (20)  do  not  change.  Therefore  they  only  have  to  be  computed  once 


and  do  not  have  to  be  revised  as  more  terms  are  added  in  series  approximation 
(13) ,  i .e. ,  larger  N. 

FQUALITY  OF  PROBABILITY  DENSITY  FUNCTION  MOMENTS 

A  very  important  property  of  expansion  (13)  is  obtained  as  follows: 


b  Q  (u) 
n  n 


(22) 


for  0  <  k  <  N 


where  we  used,  in  order,  (13),  the  orthonormality  of  (19),  and  (20).  But 
since  is  a  k-th  order  polynomial,  relation  (22)  states  that  approximation 
PN  has  exactly  the  same  moments  as  given  probability  density  function  p, 
from  order  0  through  order  N.  This  matching  of  moments  between  probability 

density  functions  p^  and  p  has  been  achieved  regardless  of  the  weighting  w 
and  its  particular  parameter  values.  Furthermore,  (22)  holds  independently  of 
whether  the  weighting-moment  equalities  in  (12)  are  satisfied  or  not. 
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The  cumulative  distribution  function  corresponding  to  approximation  pf 
is  defined  as 


u. 


N 


u 


PN(u)  =  dt  pN(t)  = 


n=0 


dt  w(t)  Qn(t)  . 


(23) 


Its  utility  depends  on  getting  closed  forms  and  simple  recursions  for  the 
general  integral  on  the  right-hand  side. 

PARAMETERS  OF  GIVEN  PROBABILITY  DENSITY  FUNCTION  p 


The  moments  of  p  were  defined  in  (3).  It  is  useful  to  define  three 
important  parameters  of  p: 

Area  A  =  j*du  p(u)  =  Pq  (p^  >  0,  but  need  not  be  1); 


Mean  Location  M  = 


_  [du  u  p(u)  _ 

"  Jdu  p(u)  "  uQ 


RMS  Width  R  = 


r  ? 

jdu(u-M)  p(u) 

(du  p(u) 


J0  VO 


(24) 


7  7 

(Conversely,  pq  =  A,  p-^  =  AM,  p£  =  A(M  +R  ).)  These  parameters 
depend  on  the  probability  density  function  p  that  we  are  trying  to  approximate 
and  can  be  computed  from  the  available  information  (9).  They  are  useful  for 
determining  where  the  major  concentration  of  p(u)  lies  on  the  u-scale,  and 
have  obvious  physical  interpretations. 


In  terms  of  the  cumulants  of  p  defined,  in  (5)-(8),  we  have  the 
alternative  expressions 
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'/a 


A  =  expQ^),  M  =  XL,  R  =^2  > 


(25) 


or  conversely 


?W"*o=*A’  ^  =  ^-©-r2- 


(26) 


GENERAL  RESULTS  FOR  THREE  LOWEST-ORDER  POLYNOMIALS  Qn 

The  weighting  function  w  and  associated  orthonormal  polynomials  satisfy 
the  following  equation: 

Jdu  w(u)  Qk(u)  Qn(u)  =  6kn  .  (27) 


Also  we  have  weighting  moments 

=  j'du  un  w(u),  with  vq  =  1  . 

It  is  then  a  straightforward  matter  to  evaluate  the  three  lowest-order 
orthonormal  polynomials: 


(28) 


Qq(u) 


Qi(u) 


Q2(u) 


=  1  , 

=  u— V i)  , 


2  L 


2/  2x 

u  (v2-Vl) 


-  u(v3-v2v1)  +  (vrv2}]  ’ 


(29) 


where 
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(30) 


The  general  expansion  coefficients  in  (20)  then  become 


(31) 


All  these  results  above  are  general  and  make  no  presumption  about  weighting 
moment  equalities  such  as  (12). 

SPECIAL  CHOICES  OF  WEIGHTING  PARAMETERS 

Suppose  that  weighting  w  has  free  parameters  that  can  be  varied  so  as  to 
make  the  mean  locations  of  w  and  p  coincide  (see  (24));  that  is. 


( 32A) 


let 


(The  reason  for  the  discrepancy  with  (12)  is  that  we  have  set  Vq  =  1  but 
have  allowed  v»q  ^  1. )  Inspection  of  (31)  gives  the  following: 


then  b^  =  0  and 


(32B) 
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Conversely,  (31)  shows  that  requiring  b-^  =  0  forces  the  choice  in  (32A)  for 
V|.  Thus  equality  of  the  first  weighting  moment  of  w  with  the  first 
(normalized)  moment  of  probability  density  function  p  implies  (and  is  implied 
by)  the  vanishing  of  the  first  expansion  coefficient  b^.  This  may  or  may 
not  be  a  useful  choice,  but,  whether  adopted  or  not,  has  no  bearing  on  the 
equality  of  probability  density  function  moments  already  demonstrated  in  (22). 


As  a  second  special  choice,  suppose  that  weighting  w  has  enough  free 
parameters  that  we  can  vary,  so  as  to  make  the  mean  locations  and  rms  widths 
of  w  and  p  coincide  (see  (24));  that  is 

u2 


P1 

let  v,  =  —  and 
1  u0 


v9  = 

2  u, 


0 


( 33A) 


(Again  we  have  used  \>q  =  1. )  Manipulation  of  (31)  yields  the  following 
conclusion: 

then  b-^  =  0  and  bp  =  0  .  (33B) 


Conversely,  imposition  of  (33B)  implies  the  results  in  (33A),  as  may  be  seen 
by  reference  to  (31).  (The  apparent  additional  solution 
\>2  =  y,/u Q  =  would  yield  an  impulse  for  w  and  is  disallowed.)  Thus 
equality  of  the  first  two  weighting  moments  of  w  with  the  first  two 
(normalized)  moments  of  probability  density  function  p  implies  (and  is  implied 
by)  the  vanishing  of  the  first  two  expansion  coefficients  b^  and  bp.  This 
common  choice  of  weighting  parameter  values  can  be  made  if  desired,  but  is  not 
necessary  (or  recommended)  for  series  approximations  by  orthonormal 
polynomials.  The  equality  of  probability  density  function  moments  in  (22) 
will  hold  whether  (33)  is  true  or  not. 
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EXAMPLE  OF  DIVERGENT  ERROR  INTEGRAL  FOR  b  =  0,  b2  =  0 


As  a  demonstration  of  what  forcing  expansion  coefficients  b^  and  b2 
equal  to  zero  can  do,  consider  probability  density  function 


2,  2n 


/  \  2  u  exp(-u  / u  )  ^ 

P(u)  -  f°ru> 


( y  >  -1»  to  >  0) 


with  moments 


(34) 


(35) 


This  class  of  probability  density  functions  includes  the  one-sided  Gaussian, 
Rayleigh,  and  Maxwell  as  special  cases,  for  y  -  0,  1,  2,  respectively. 


Consider  also  the  Hermite  weighting  given  in  (1),  which  has  moments 
(11)  equal  to 


Vq  =1,  =  a,  \>2  =  .  (36) 

If  we  now  insist  on  property  (33B),  then  (33A)  yields 


But  the  leading  integral  in  minimum  error  E^  in  (21)  is  convergent  only  if 

p 

p  (u)/w(u)  decays  sufficiently  rapid  for  large  u.  We  have  from  (34)  and 
(1),  the  dominant  behavior 


p  (u)/w(u)oC  expf-  -r—  +  —0  \  for  large  positive  u  , 


(38) 


wherecC  denotes  proportionality,  but  disregards  the  exact  scale  factor  and 
subdominant  behavior.  Thus  the  integral  in  (21)  is  convergent  only  if 
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U 


2(y+1) 


(39) 


However,  calculation  of  (39)  reveals  that  this  inequality  is  never  satisfied 
for  any  value  of  y  >  -1,  the  function  on  the  right-hand  side  starts  at  0  when 
y  =  -1,  and  increases  monotonical ly  towards  1  as  y  >  +o0>  behaving  like 
1  -  1/ ( 4y)  in  this  limit. 


Thus  expansion  of  probability  density  function  (34)  according  to  a 
Hermite  weighting  has  an  infinite  error  integral  (21)  (and  perhaps  a  divergent 
series  expansion)  regardless  of  the  values  of  y  and  u>  in  the  true  probability 
density  function,  if  we  insist  on  expansion  coefficients  b^  =  b2  =  0.  Yet 
if  we  relax  requirement  (33B),  and  choose  3  according  to  (39)  such  that 
B  >  u/2,  the  error  integral  in  (21)  is  certainly  finite,  regardless  of  a. 


However,  making  the  error  integral  in  (21)  finite  is  not  the  whole  story, 
in  so  far  as  realizing  useful  approximations  to  the  probability  density 
function  or  cumulative  distribution  function.  An  example  of  probability 
density  function  (34)  was  taken  with  y  =  3,  u  =  1.  When  a  and  b  were  chosen 
according  to  (33)  and  (37)  (giving  a  =  .48  <  .5  =  u/2),  the  expansion 
coefficients  jb 3  initially  decreased  to  approximately  IE— 3  at  n  =  40  terms, 
and  then  diverged;  yet  a  plot  of  the  approximate  exceedance  distribution 
function  obtained  by  a  Hermite  expansion  overlaid  the  exact  answer  down  to  the 
IE— 16  level.  On  the  other  hand,  when  the  weighting  parameters  in  the  Hermite 
expansion  were  chosen  as*  a  =  0,  b  =  .7  >  .5  =  u/2,  giving  b^  ^  o  and 

*This  is  example  B  in  a  later  section 
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b2  4  0,  the  expansion  coefficients  {b  j  decreased  to  the  IE— 4  level  at 
n  =  70  before  round-off  error  became  dominant;  despite  this  apparent 
improvement  in  coefficient  level,  the  approximate  exceedance  distribution 
function  overlaid  a  plot  of  the  exact  result  down  to  the  IE-10  probability 
level,  which  is  several  orders  of  magnitude  worse  than  the  above  result.  Thus 
emphasis  on  getting  a  convergent  error  integral  in  (21)  may  not  always  be 
desired. 

For  Hermite  weighting  (1)  and  the  class  of  probability  density  functions 
which  decay  as  exp(-uq)  as  u  >  +<*>,  the  error  integral  is  always  convergent 
if  q  >  2,  and  always  divergent  if  q  <  2.  So  an  exponential  probability 
density  function,  like  uY  exp(-u/u>)  for  u  >  0,  always  yields  a  divergent 
error  integral  when  expanded  in  a  Hermite  series. 

For  generalized  Laguerre  weighting  (2),  it  is  necessary  to  consider 
u  =  0+  and  u  =  +«>  separately.  If  probability  density  function  p  behaves  like 
uY  as  u  >  0+,  then  a  finite  error  integral  requires  that  we  choose 

a  <  1  +  2y.  Coupled  with  the  finite  area  restriction  on  weighting  w,  a  range 

of  values  of  a  is  allowed,  namely,  -1  <  a  <  1  +  2y,  this  range  always  exists 

since  y  >  -1  is  necessary  for  the  probability  density  function  itself  to  have 

finite  area. 

If  also  the  probability  density  function  behaves  as  exp(-u/w)  as  u  >  +<»  , 
then  a  finite  error  integral  with  generalized  Laguerre  weighting  requires  that 
we  choose  a  >  u>/ 2.  So  the  range  of  choice  of  a  is  open  on  the  large  side, 
whereas  that  for  a  is  a  limited  one,  for  this  particular  class  of  probability 
density  functions. 
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HERMITE  EXPANSION 


In  this  section,  we  will  deal  exclusively  with  weighting  (1), 
w(u)  =  j  <t>  for  all  u 

where 

“A  2 

x )  =  (2ir)  exp(-x^/2). 

This  weighting  has  two  free  parameters,  a  and  a,  and  moments 

Vg  =1,  vi  =  v2  =  ’  (42) 

If  and  ^2  are  specified,  the  parameters  must  then  satisfy  a  =  v^, 

2  1/2 

B  =  (v2  -  vi)  '  •  However,  we  shall  keep  a  and  b  general  and  unspecified. 

PROPERTIES  OF  POLYNOMIALS  AND  EXPANSIONS 


(6  >  0)  , 


(40) 


$(x)  =  j  dt  d(t)  . 


(41) 


The  orthononnal  polynomials  associated  with  weighting  (40)  are  the 
Hermite  polynomials  [5,  22.1.2  and  22.2.15] 


Q  (u)  =  He  (nl) 

n  n  \  b  / 


-A 


for  n  >  0 


(43) 


The  expansion  coefficients  are  given  by  (20)  as 

t>n  =  j"  du  p ( u )  Qn(u)  =  (nL)-1^2  cp  for  n  >  0  , 


(44) 
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where  we  define 


:n  -  J  du  p(u)  Hen(iyL)  for  n  >  0  . 


(45) 


The  approximate  probability  density  function  then  follows  from  (13)  in  the  form 

N  N 

PN(u)  =  w(u)  ^  bn  Qn(u)  =  I  an  ^nfnr)  ’  (46) 


n  n 


6  \3 


n=0  n=0 

where  we  used  (40),  (43),  (44),  and  defined 


,  , xl/2  ,  ,  , x-1/2 

(n')  an  =  bn  =  (ni) 


n  n 


c  for  n  >  0  . 
n  — 


(47) 


These  three  different  coefficients  in  (44)-(47)  are  introduced  for  convenience 
in  further  equation  manipulations.  Expansion  coefficient  bp  is  the 
geometric  mean  of  auxiliary  coefficients  an  and  cn  (with  polarity). 

Expansion  (46)  is  also  called  a  Gram-Charl ier  series  of  type  A  [2,  p.  222], 

[3,  p.  156]. 


The  approximate  cumulative  distribution  function  corresponding  to  (46)  is 

PN(U)  =  f  it  PN(t)  =  r  *  (if1)  Hen  {¥)  ■ 


n=0 


T 


n=0  — <*> 


n  f  dx  <*(x)  Hen(x)  =  ao^(T)  -  «S(T)  an  He^T)  ,  (48) 

-oo  n=l 


where 


and  we  used  (41)  and  [5,  22.11.8]. 


(49) 
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The  Hermite  polynomials  {He^  satisfy  the  recurrence  [5,  22.7.14] 


Hen(x)  =  x  Hen_^(x)  -  (n-1)  Hen_2(x)  for  n  >  2  ,  (50) 


with  starting  values  He0(x)  =  1,  He-^x)  =  x  [5,  22.3.11].  The  highest- 
order  term  in  Hen(x)  is  x11,  with  coefficient  1  [5,  22.1.2  and  22.3.11]. 

The  magnitude  of  the  term  multiplying  bn  in  (46)  has  an  envelope  that  decays 
approximately  as  n-^4  with  n,  regardless  of  argument  u.  This  may  be  seen 
by  using  (47)  and  (49)  to  get 


Hen(T)  =  bn(ni) 


-1/2 


Hen(T)oC  b 


n(n"+1'2  e-")'1'2  (n/e) 


n/2 


=  bn  n 


-1/4 


as  n  >+<»,  for  all  T,  (50A) 

where  we  also  used  [5,  6.1.39  and  22.5.18]  and  [7,  8.22.8].  Here,oC  denotes 
proportionality  and  we  have  taken  the  magnitude  of  the  terms;  the  exact  scale 
factor  of  proportionality  will  be  presented  in  a  later  section  where  the 
errors  of  the  approximations  are  estimated.  So  if  bn  were  to  decay  faster 

O  /  A 

than  n  1  ,  the  probability  density  function  series  in  (46)  would  converge 
absolutely. 


Conditions  are  better  for  the  cumulative  distribution  function  series  in 
(48);  namely,  based  on  the  above  result,  there  follows  (for  the  envelope) 

an  Hen-1(T)  =  bn(nl)"1/2  Hen-1(T)  =  bn  C(n-l):3“1/2  He^T)  = 

c C  bn  n“^  n-^4  =  bn  n-^^4  as  n  >  +0°,  for  all  T  .  (50B) 

Thus  if  b^  decays  faster  than  n-^4,  the  cumulative  distribution  function 
series  converges  absolutely.  Furthermore,  if  the  leading  error  integral  in 
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o 

(21)  is  finite,  the  sum  of  bp  must  be  finite,  meaning  that  bn  must  decay 
faster  than  n-^.  So  we  can  conclude  that  if  the  error  integral  is  finite, 
the  Hermite  series  for  the  cumulative  distribution  function  in  (48) 
converges.  (Notice  that  this  particular  decay  n-^  0f  bn  is  not 
sufficiently  fast  to  make  the  same  conclusion  about  the  Hermite  series  for  the 
probability  density  function  in  (46).)  The  above  are  sufficient  conditions  on 
expansion  coefficients  {b^,  and  are  not  necessary. 

EXPANSION  OF  CHARACTERISTIC  FUNCTION  f 


The  coefficients  an  and  cn  were  defined  in  (45)  and  (47).  Then  the  sum 
2  an  w"  “  ST  cn  w"  -j>  nT  Idu  »(u)  Hen  (tO  = 


n=0 


n=0 

CO 


n=0 


.  J  du  p(u)  He„  =  jdu  p(u)  exp^jS  „  -  j 


n=0 

=  exP  (-  \  W^  -f  w\  f  (fj  , 


(51) 


where  f  is  the  characteristic  function,  and  where  we  used  (45),  [5,  22.5.19 

and  22.9.17],  and  (5).  Letting  w  =  ez,  we  have 

00 

(52) 

n .  u 

n=0  n=0 


(  \  °° 

f  (z)  exp  \-az  -  j  e2z2J  =  an  (ez)n  =  cn  (6z)n 


Thus  {a^  and  {cp^j  are  the  coefficients  in  these  power  series  expansions 
of  the  function  f(z)  exp(-az  -  b^z^/2),  where  f  is  the  characteristic 
function  corresponding  to  probability  density  function  p,  and  a  and  e  are 
arbitrary.  A  special  case  of  (52)  is  given  in  [2,  17.6.10]. 
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Collecting  (46)  and  (52)  together  for  comparison,  and  assuming  that 
Pn  >  p  as  N  >  +c*>,  we  have 


(53) 


Thus  expansion  of  probability  density  function  p  in  an  infinite  Hermite  series 
is  equivalent  to  an  expansion  of  a  modified  form  of  the  characteristic 
function  in  a  power  series,  according  to  (53).  Equations  ( 51) — ( 53)  will  serve 
as  very  convenient  starting  points  for  the  derivation  of  several  alternative 
recurrences  for  the  expansion  coefficients  Notice  that  weighting 

parameters  a  and  a  are  completely  unrestricted  in  (52)  and  (53),  except  that 

3  >  0. 

An  analogous  result  holds  for  N  finite,  but  must  be  derived  in  a 
different  fashion,  because  we  no  longer  can  use  infinite  sum  [5,  22.9.17]. 
Define  the  Fourier  transform  of  (46)  as  the  N-th  order  approximation  to  the 
characteristic  function: 


N 


n=0 


N 


n=0 


N 
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N  f 

exp(ia^)  an  (i3?)n  ]  dt  exp(iest)  d(t)  = 

n=0 


exp  ^iaf  +  \  e2(i?)^) 


N 


an  (ie?)n  , 


(54) 


n=0 


where  we  used  [5,  22.11.8]  in  line  4,  and  repeated  integration  by  parts  in 
line  5.  This  result  is  the  leading  N  terms  of  (53).  As  a  by-product  of  this 
derivation,  we  have 


j~dt  exp(zt)  d(t)  Hen(t)  =  exp^|  z2^  zn  . 


(55) 


COEFFICIENTS  RECURSIVELY  VIA  CUMULANTS 


We  are  now  in  a  position  to  obtain  some  useful  recursive  relations  for 
the  expansion  coefficients  {a^  in  ( 51) -( 54) .  The  first  one  is  obtained  by 
taking  the  JU  of  (51) : 

Jw  f  (j)  -  "I  w  -  7  w2  -  JtM  an  A  (56) 

(n=0  J 


Then  using  (7)  and  identifying  the  right-hand  side  of  (56)  as  a  new  power 
series,  we  have 
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There  follows  immediately 

(% 


h  =< 
n  \ 


n!  6 


B 


1 

2R 

B 


for  n  4  1,2 


for  n  =  1 


(58) 


-  1  for  n  =  2 


v_ 


But  equality  of  the  right-hand  sides  of  (56)  and  (57)  also  requires  that 

c*>  r  sl  "''i 


a  wn  =  exp 
n 


n=0 


(i 

[n=0 


(59) 


It  is  shown  in  appendix  A  that  a  recursive  solution  to  (59)  for  the  ^a  ^  is 
given  by 

n 


“n  =  n  2  m  hm  an-ra  for  "  -i  aO  *  exP<hO> 


(60) 


m=l 


Then  eliminating  ^h^  by  means  of  (58), 

(?  -  f)  3n-l  +  (fz  ~  an-2  +  2 


an  =  n 


X, 


(m-i) ;  e 


m  n-m 


for  n  >  1, 
aQ  =  expC^),  (61) 


where  an  =  0  for  n  <  0,  and  the  sum  is  zero  for  n  <  3. 


Now  define  normalized  cumulants  (excluding  n=0)  according  to 


-y  X, 
n  (n-1):  s" 


for  n  >  1 


(62) 
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Then  (61)  becomes 


n 


a 


_  1 
n  n 


m=3 


aQ  =  expOtg)  .  (63) 


This  convolution  is  the  desired  recursion  for  expansion  coefficients  ^an] 
via  cumulants. 

As  particular  cases,  we  have 


( 64A) 


Parameters  a  and  6  (>0)  are  completely  arbitrary  in  the  above  three  equations, 


available  cumulants  of  the  probability  density  function 


under  consideration. 


Observe  that  if  we  choose  a  =  X\  =  M  and  b  =  R  (see  (24)-(25)), 


which  is  a  very  common  choice,  we  have  a  =  0  and  ^  =  0;  this  is  a 
special  case  of  the  general  property  (33)  stated  earlier.  This  special  choice 
of  a  and  b  corresponds  to  choosing  the  mean  location  and  rms  width  of  Hermite 
weighting  (40)  identical  to  those  same  parameters  of  the  given  probability 
density  function.  There  then  also  follows,  in  this  special  case, 


for  n  >_  6  when  ct  =  ^ >  B  =  •  (64B) 
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COEFFICIENTS  DIRECTLY  VIA  MOMENTS 


Before  we  begin  this  derivation,  we  present  the  following  useful 
expansion  [5,  22.9.17  and  22.5.19]: 

\  oo 

exp  (-  \  y2  +  xy)  =  Hen(x)  yn  . 


(65) 


n=0 


We  now  again  refer  to  (51)  and  expand  the  terms  as  follows: 


n=0 


oC 


um  (? 


m 


bHe k(-f)"k  2  ST 

k=0  m=0 

where  we  utilized  (65)  and  (6).  Equating  coefficients  of  wn  on  both  sides 
of  this  equation,  we  have 


(66) 


an  " 


n-k 


k=0 


(n-k) I  b 


n-k 


for  n  >  0 


(67) 


We  now  define,  for  convenience,  the  normalized  Hermite  polynomials 


He  (x)  =  —7  He  (x)  for  n  >  0  ,  (68) 

n'  '  n.‘  n  '  —  ’ 

and  the  normalized  moments 

C  =  — - —  for  n  >  0  .  (69) 

n  ,  „n  — 

nl  b 

(Notice  the  difference  with  the  definition  of  the  normalized  cumulants  (62).) 
Then  (67)  becomes 
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a 


n 


for  n  >  0  , 


(70) 


which  gives  expansion  coefficients  [an]  directly  in  terms  of  the 
(normalized)  moments  of  the  given  probability  density  function.  The 
recurrence  in  (50)  can  be  used  to  generate  the  Hermite  factors  needed  in 
convolution  (70).  Parameters  a  and  6  (>0)  of  weighting  (40)  are  arbitrary. 
Heo(x>  l,  Hen(x)=  ^  [x  Hen_a(*)]  -for  a^2. 

As  particular  cases,  we  have 


a0  -  Uq’ 


al 


V1  ~  ap0 

6 


2aU]  +  (ot^-0^)  ug 


2fi 


(71) 


These  agree  with  (64)  which  utilized  curnulants.  If  we  make  the  special  choice 

2  2 

of  a  =  o-j/pq  and  0  =  v*2/u0  “  ’  tden  ai  =  ®  and  a2  = 


An  alternative  more  direct  derivation  of  (67)  is  possible,  from  (47), 
(45),  (B-3)  in  appendix  B,  and  (3), 

an  -  h  cn  '  TTT  [  du  P(u>  Hen(ir)  = 

.  f"  )  He„  f-  k 

k=0 


2-  Hek(‘  s)  („-k):ke"-k 


for  n  >  0  . 


(72) 
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COEFFICIENTS  RECURSIVELY  VIA  MOMENTS 


Before  we  begin  this  derivation,  we  replace  x  >  -ix,  y  >  iy  ■ 

exp(|  y2  +  xy)  -  jf  ±t  "'n1'")  (iy>n  =  iT  Hin(x>  7 

n=0  n=0 

where  Hin(x)  is  a  real  n-th  order  modified  Hermite  polynomial  in  : 

Hi n ( x }  =  in  Hen(-ix)  for  n  2  0  • 

The  recursion  for  these  polynomials  follows  immediately  from  (50) 
Hin(x)  =  x  Hin_1(x)  +  (n-1)  Hin_2(x)  for  n  2  2  , 


with 


starting  values  HiQ(x)  =  1,  Hi^x)  =  x.  The  difference  with 

2 

the  polarity  of  the  last  term;  thus  for  example,  Hi2(x)  =  x  +  1, 
Hi2(x)  =  x°  +  3x,  versus  HegU)  =  x^  -  1,  He-^x)  =  x^  -  3x. 


We  now  rewrite  (51)  in  the  following  form: 

f(?)  -  exp(?  1,2  +  f  ")  J.  ara  w"  ' 

m=0 

Expanding  in  power  series  by  means  of  (6)  and  (73), 


<*> 


n=0 


nT  un (&)  ’  S  iC-H1k(r  "k^. 

k=0  m=0 


a  w 
m 


Equating  coefficients  of  w  ,  there  follows 


n!  e 


k=0 


lrHik(f)  an-k  for  n2°  , 


or 


in  (65)  to  get 
1  ,  (73) 

<  defined  by 

(74) 
as 

(75) 

(50)  is 

(76) 

(77) 

(78) 
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A 


for  n  >  0  , 


where  we  have  used  normalized  moments  (69),  and  defined  the  normalized 
modified  Hermite  polynomials 

Hin(x)  =  'pjy  Hin(x)  for  n  >  0  . 


(79) 


(80) 


Finally,  the  desired  recursion  for  expansion  coefficients  ja^  in  terms  of 
the  moments  follows  as 


A 


n 


for  n  >  0  . 


Parameters  a  and  0  (>0)  are  arbitrary  in  (81)  and  (69). 


(81) 


SUMMARY 


The  approximations  to  the  probability  density  function  and  cumulative 
distribution  function  are  given  by  (46)  and  (48),  respectively,  where  a  and  & 
are  arbitrary  constants,  except  that  &  >  0.  The  functions  <t>  and  $  are  defined 
in  (41),  while  the  Hermite  polynomials  (He^  are  available  via  (50).  The 
expansion  coefficients  [a^  are  given  by  the  three  alternatives  (63),  (70), 
(81),  in  terms  of  normalized  cumulants  (62),  normalized  moments  (69), 
normalized  Hermite  polynomials  (68),  and  normalized  modified  Hermite 
polynomials  (80)  and  (74).  Programs  for  all  three  alternative  procedures  for 
determining  expansion  coefficients  [an]  are  listed  in  an  appendix.  The 
basis  for  these  relations  is  the  characteristics  function  expansion  in 
(51) -( 53) . 
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GENERALIZED  LAGUERRE  EXPANSION 


This  section  will  treat  weighting  (2),  namely. 


for  u>0  (a  >  -1,  e  >  0) . 


(82) 


This  weighting  is  a  special  case  of  the  three-parameter  weighting 


(u-y)“  exp 


for  u  >  y  , 


(83) 


sa+1  run) 


which  is  the  most  general  scaled  linear  shift  of  the  generalized  Laguerre 
weighting  [5,  22.2.12] 


xa  exp(-x)  for  x  >  0  . 


(84) 


We  will  consider  only  y  =  0  here.  For  a  probability  density  function  Pq(u) 
which  is  known  to  be  nonzero  only  for  u  >  Ug,  we  would  consider  the  modified 
probability  density  function  p(u)  =  p0(u+ug),  because  then 
p(u)  4  0  only  for  u  >  0,  and  the  simpler  weighting  (82)  would  be  directly 
applicable.  This  procedure  is  equivalent  to  choosing  y  =  uQ  in  the 
three-parameter  weighting  (83)  above,  and  requires  knowledge  of  Ug.  We 
presume  that  p(u)  ^  0  only  for  u  >  0  henceforth  in  this  section,  and  that  any 
necessary  shifting  has  already  taken  place. 

Weighting  (82)  has  two  free  parameters,  a  and  e,  and  moments 


v  =  ( a+l)  Bn  for  n  >  0  . 
n  '  'n  — 


(85) 
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In  particular. 


(86) 


If  v-^  and  \>2  are  specified,  then  the  parameters  must  satisfy 


2 


a 


- 2  -1  .  e 

v2  ”  V1 


(87) 


However,  we  shall  keep  a  and  e  general  and  unspecified  except  for  the 
conditions  in  (82). 

PROPERTIES  OF  POLYNOMIALS  AND  EXPANSIONS 

The  orthonormal  polynomials  associated  with  weighting  (82)  are  the 
generalized  Laguerre  polynomials  [5,  22.1.2  and  22.2.12] 


for  n  >  0,  u  >  0  . 


(88) 


The  expansion  coefficients  are  given  by  (20)  as 


for  n  >  0  , 


(89) 


where  we  define 


for  n  >  0  . 


(90) 


The  approximate  probability  density  function  follows  from  (13)  according  to 
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N 


pN(u)  =  wiu)  J>  b  Q  (u)  = 

IN  n=0  n 

_  ua  exp(-u/ b)  ^  L^aY- 

'7Trun)  u  n  nU' 


where  we  used  (82),  (88),  (89),  and  defined 

v  1 

n  l  \ 


ni 


n( 


for  u  >  0  , 


(91) 


c„  for  n  >  0  . 


(92) 


These  three  different  coefficients  in  (89)-(92)  are  introduced  for  convenience 
in  further  equation  manipulations.  Expansion  coefficient  b^  is  the 
geometric  mean  of  auxiliary  coefficients  an  and  cn  (with  polarity). 


The  approximate  cumulative  distribution  function  corresponding  to  (91)  is 

u  u 

r  N 

PN(u)  =•  j  dt  pN(t)  =  a 


o 


N 


f  ,+  ta  exp(  -t/g)  .  (a)ft\  _ 


i 

r(  «+d 


n=0 


an  !n (?)  for  u  >  0  ’ 


(93) 


where  we  define 


I^(y)  =  ^dx  xa  e  x  lJ“)(x)  for  n  _>  0,  y  >  0  . 


(94) 


These  quantities  are  evaluated  in  appendix  C;  when  substituted  in  (93),  they 
yield  (95) 


v  a+l 


p  (u/B)  exp(-u/B) 

Vuj  “  r(T+l) 


N  a 


0  F  ( l  •  ^  u  1  -*■  ^  n  |  ( ®+l) AA 

a+T  lFlU’a 2V  n  L  n-1  U/ 


for  u  >  0, 


where  jF ^  is  the  confluent  hypergeometric  function. 
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The  generalized  Laguerre  polynomials  satisfy  the  recurrence 

[5,  22.7.12] 

L^(x)  =  ^  j(a-l+2n-x)  L^“|(x)  -  (a-l+n)  L^(x7J  for  n  >  2  ,  (96) 

with  starting  values  L^(x)  =  1,  L'“'(x)  =  a+l-x  [5,  22.4.7].  The  highest 
order  term  in  L  ^  (x)  is  ( — x ) n / n I  [5,  22.1.2  and  22.3.9];  this  is  distinctly 


different  from  the  coefficient  1  for  the  Hermite  polynomials.  Yet  the 
envelope  decay  with  n  of  the  generalized  Laguerre  series  for  the  probability 
density  function  and  cumulative  distribution  function  are  identical  to  those 
of  the  Hermite  series,  for  u  >  0.  To  prove  this,  use  (91)  and  (92)  to  get 


as  n  >  +<*>,  for  u  >  0  ,  (97) 

where  we  also  used  [5,  6.1.39]  and  [7,  8.22.1].  Again,  oCdenotes 
proportional ity  with  n  only;  the  exact  scale  factor  will  be  presented  in  a 
later  section  where  the  errors  of  the  approximations  are  estimated.  So  if 
bn  decays  faster  than  n-^4,  the  probability  density  function  series  in 
(91)  converges  absolutely. 


For  the  generalized  Laguerre  series  of  the  cumulative  distribution 
function  in  (95),  we  have,  for  the  envelope  of  the  general  term, 
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1 

n 


a 


n 


.iW 

a+1  1 


L(a+1)/U\ 

n-1  U/ 


cCbn  i  (n~°f  (n-1)  2  4  ~  bn  n_3/4  as  n  >  +*,  for  u  >  0  .  (98) 

Thus  if  bn  decays  faster  than  n-4  ,  (95)  converges  absolutely.  And  if 
the  error  integral  (21)  is  finite,  this  property  of  the  [b^  is  true.  So  if 
error  integral  (21)  is  finite,  the  generalized  Laguerre  series  for  the 
cumulative  distribution  function  converges  absolutely  for  u  >  0;  this  is  a 
sufficient,  but  not  necessary,  condition. 


For  zero  argument,  the  generalized  Laguerre  polynomials  behave 
differently  for  large  n.  From  [5,  22.4.7  and  6.1.39], 

(«+l), 


i  ( o) / q\  /n+«\  a  n  -  na 

L  n  (0)  =  (  n  )  =  ~nT~ 


as  n  >  +°o. 


(99) 


Then  (97)  and  (98)  are  both  replaced  by  bn  na/  as  n  >  +00.  However,  for 
a  >  U,  the  probability  density  function  in  (91)  is  zero  at  u  =  0  due  to  the 
u°  term,  so  there  is  no  need  to  perform  the  sum  then.  And  the  cumulative 
distribution  function  is  always  zero  at  u  =  U,  again  eliminating  the  need  to 
evaluate  the  sum  in  (95).  So  the  difference  in  behavior  at  u  =  0  is  of  no 
consequence. 
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EXPANSION  OF  CHARACTERISTIC  FUNCTION  f 


The  coefficients  and  cn  for  the  generalized  Laguerre  series  were 
defined  in  (90)  and  (92).  Then  the  sum 


o C 


n=0 


cn  w  = 


n=0 

OO 


•T. 


w"  l  du  p ( u )  L(“}(u/0)  = 
J0 


=  f  du  p(u)  wn  L^(u/b)  =  f  du  p(u)  (1-w)  “  1  = 

Jo  n=0  n  J0  V  1-w/ 

-  a-w)— 1  f(^)  , 


(100) 


where  f  is  the  characteristic  function,  and  where  we  used  (90),  [5,  22.9.15], 
and  (5).  Thus  [cn]  are  the  expansion  coefficients  of  the  right-hand  side  of 

(100)  in  powers  of  w.  If  we  let  w  =  ’  we  l1376  the  expansion  for  the 

characteristic  function 


f(z) 


d-62) --1 


O0 


n=0 


-6Z 


n  Vl— 6 z 


(101) 


corresponding  to  given  probability  density  function  p.  Weighting  parameters  a 
and  b  are  arbitrary  in  (100)  and  (101). 


Collecting  (91)  and  (101)  together  for  comparison,  and  assuming  that 
P[vj  >  P  as  N  >  +  0°,  we  have,  upon  use  of  (92), 


P(u) 


ua  exp(-u/B )  "C  , (a)/u\ 

nTT)  U  "  "  (») 


for  u  >  0  , 


f(i?)  =  (l-iBf)"0-1 


©O 


n=0 


(“+1  )r 
FTi 


(102) 
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Thus,  expansion  of  probability  density  function  p  in  an  infinite  generalized 
Laguerre  series  is  equivalent  to  an  expansion  of  the  corresponding 
characteristic  function  in  the  series  of  the  particular  form  in  (102). 
Equations  ( 100 ) — ( 10 2)  will  serve  as  very  convenient  starting  points  for  the 
derivation  of  several  alternative  recurrences  for  expansion  coefficients 
^a^  „  We  reiterate  that  a  and  &  are  arbitrary  in  the  above,  except  that 

a  >  -1,  6  >  0. 


An  analogous  result  holds  for  N  finite,  but  must  be  derived  differently 
since  we  can  no  longer  use  infinite  sum  [5,  22.9.15].  Define  the  Fourier 
transform  of  (91)  as  the  N-th  order  approximation  to  the  characteristic 
function : 

fN(if)  =  j  du  exp(ifu)  pN(u)  = 


00 

f  j  / \  u  exp( -u/b) 
=  du  exp(ifu)  -prtr - - 

o  run) 


N 

S 

n=0 


a  L(  “Y- 
an  L  n  U 


fUn  2  an  fdt  expll6ft)  e_t  Lln)(t) 


(103) 


In  appendix  0,  it  is  shown  that 

iwt  o  -t  ,(o)m  P(  a+  l+n )  ( — i  oj  ) n 

dt  e  t  e  L  (t)“  Pi  TT~Unn 


(104) 


Substitution  in  (103)  then  yields 


fN(i5) 


a  (“n)n  (-ief)n 

"  n;  ( 1-i 6? ) a+l+n 


1 -leri" 

(l-16|)“n+n  ’ 


(105) 


where  the  last  relation  follows  by  use  of  (92).  This  result  is  the  leading  N 
terms  of  (102). 
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COEFFICIENTS  RECURSIVELY  VIA  CUMULANTS 


We  can  now  obtain  some  useful  recursive  relations  for  expansion 
coefficients  and/or  in  (100)-(105).  We  start  by  taking  the 
of  (100): 


J?n  l  ^  c 
ln=0 


n  wnl  =  -(a+1)  in(l-w)  +  J>n  f  • 


(106) 


Identify  the  left-hand  side  as  a  new  power  series,  and  use  (7)  and  [5,  15.1.8] 
to  yield 

,  k 


Ihnw".(aH)  ill,"*  2 

n=0  n  n=l  n  k=0  K*  K  V1  w  . 


/  ^  1  n  ^  ^  (  A  k  S  k  m  m 

-  (“+1)  2i  +  <>-  - r—  w  —r-  w  . 

n=l  n  k=0  kl  bk  m=0  m‘ 

Equating  coefficients  of  wn,  there  follows  hg  =  %g,  while  for  n  _>  1, 

„  1,  ...  +  A  (k)n-k 

n  "  "  “  m  i:  Bk  (n-k):  ' 


■  i  h*  ±  <-»k  (!)* 


(107) 


(108) 


where  we  used  the  normalized  cumulants  defined  in  (62). 


But  since  the  left-hand  sides  of  (106)  and  (107)  are  equal,  we  have 


n  W"  =  6XP  I  n  h"  Wf  ’ 

n=0  Ln=0  J 


(109) 


or  via  appendix  A,  the  recurrence 
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cn  n  m=I  111  Ch-<" 


for  n  >  1,  c  =  exp(h0 ) 


(110) 


Finally,  define 


d  =  m  hK  for  m  >  1 

m  *»  — 


(111) 


far  notational  convenience  and  thereby  obtain 


m 


dra-“n+|g  forl"i1' 


1  11 

c  =  —  d  c 
n  n  m  n-m 


for  n  >  1,  c0  =  exp(^)  ,  (112) 


by  means  of  (108)  and  (110),  respectively.  Equation  (112)  is  a  recursive 
relation  for  expansion  coefficients  £c^  in  terms  of  cumulants  and 

auxiliary  variables  {dm] .  The  [an|  are  immediately  available  via  (92). 


As  particular  cases,  we  have,  employing  (62), 


c 


(«+2)(«+l)  -  2(a+2)^ 

P 


+ 


(113) 


Parameters  a  and  0  are  completely  arbitrary  in  all  the  above  equations,  except 
that  a  >  -1  and  0  >  0,  and  are  the  available  cumulants. 


Observe  that  if  we 


let  a+l  =  y— 
*2 


2 

U1 


y2y(Tul 


7 


and 


0  = 


*1 


1^0  ~  » 


1  R 


(114) 
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then  c^  =  0  and  c2  =  0  (here  we  also  used  (8)  and  (25));  this  is  a  special 
case  of  general  property  (33)  stated  earlier.  Since  the  probability  density 
function  p(u)  is  nonzero  only  for  u  >  0,  then  >  0  and  >  Q,  giving 
allowable  solutions  to  (114)  in  all  cases.  There  then  also  follows,  along 
with  d^  =  d£  =  0  in  this  special  case,  the  explicit  results 

cQ  =  expO^),  c1  =  0,  c2  =  0, 


+  300^  X*2  Xf  -  3QX5Z2tf  +%6  c0. 


+  42OX3  Xl7-\+  4200>^  X\  +  70)^2  X\  -  35^  ^  - 

-  630 +  42^  \X\-Y-j  cQ. 


(115) 
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These  relations  have  been  confirmed  by  numerical  comparison  with  (112). 

These  results  greatly  extend  those  of  [1,  ( 129) — ( 131) 2,  where  the 
equivalent  of  our  c^  is  given  (in  terms  of  moments  instead  of  cumulants,  and 
with  Kq  =  0),  along  with  the  comment  that  "the  higher-order  coefficients  are 
so  complicated  that  the  whole  value  of  this  type  of  series  seems  to  depend  on 
the  fact  that  the  first  term  alone  (cQ)  is  often  a  good  approximation."  We 
find,  on  the  other  hand,  that  not  only  can  we  avoid  the  special  choice  in 
(114)  and  the  corresponding  complicated  special  results  in  (115),  but  we  can 
handle  any  a,e  pair  and  get  very  high-order  coefficients  cn,  simply  by  using 
the  recurrence  in  (112),  which  is  easily  programmed.  The  only  thing  we  lose 
are  explicit  results  of  the  type  given  in  (115j;  however,  the  latter  are  so 
complicated  that  they  are  of  limited  utility  anyway. 

COEFFICIENTS  DIRECTLY  VIA  MOMENTS 


We  will  need  the  following  expression  [5,  22.3.9]: 


ni  i  (a),  x  n! 
T^rr  L  n  (x)  "  T^+T) 


n+o\  (~x)K  _  (-x) 

n-k/  k T~ 


n  k=0  **  k=0  'kl  (a+ljk  ’ 


(116) 


Then  (92),  (90),  and  (3)  yield,  for  n  >  0, 

a  =  /n‘ -t  t  c  =  ,n’ ■<  t  T  du  p(u)  L^aV— )  = 

n  lo+l)n  n  (o+l)n  n  U ) 


“•  'E)^  fduplu)  (-u/B)k  =  ±(-l)k^  “k 
k  Jo 


k=0 


k=0 


k'  (a+l)k  Bk  ‘ 


(117) 


It  is  useful,  in  this  generalized  Laguerre  series  case,  to  define  an 
alternative  set  of  normalized  moments 
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“  (a+l)  Bn 
'  n 


for  n  >  0  . 


(118) 


(Although  this  seems  to  be  very  different  from  the  earlier  normalization  in 
(69),  (118)  actually  reduces  to  (69)  for  the  a  here  equal  to  zero.)  When 
(118)  is  utilized  in  (117),  we  have  the  desired  expression  for  expansion 
coefficients  £a  } ,  directly  in  terms  of  (normalized)  moments,  in  the 
surprisingly  simple  form 

n 


(-l)k  (j)  u.  for  n  >  0  . 


an  jgjj  v  kk/  Mk 


(119) 


Parameters  a  and  3  in  (118)  are  arbitrary,  except  that  a  >  -1,  8  >  0. 


As  particular  cases,  ( 117)— ( 119)  yield 


=  A 


O’ 


U1 

al  =  y0  (a+l)B  ’  a2  "  y0 


2u- 


(a+1)6  (ct+1)  (a+2)  32 


These  agree  with  (113)  which  utilized  cumulants.  If  we  make  the  special 

2  2  2 
choices  of  a+l  =  ^/(v^Aq  -  y^)  and  3  =  (^2U0  “  ul^  (v]Uq)> 

then  a-^  =  0  and  -  0;  this  is  a  common  approach  to  the  approximation 

problem,  but  totally  unnecessary. 


(120) 


An  alternative  derivation  of  the  direct  moment  relation  (117)  is 
possible:  from  (100),  (6),  and  [5,  15.1.8], 
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n  , ,  \-a-l  v  1  /-w/b 

c  w  =  1-w  777  uu  -777 

n=0  n  k=0  k’  k  V  1_w 


06 


¥T“k  H) 


oO 


k:  “k  ('  s)  Jj 


k=0  \  p/  k=0  ^  N  ^  p/  m=0 

Equating  coefficients  of  wn,  we  have,  for  n  >  0, 


n 


1  f  nk  (a+1+l<>n-k  (<,+1)n  ■£- 


cn  kl  “kl  b  1  (n-k): 


n: 


k=0 


(-1) 


k  fn 


(a+l+k) 


m  wm 


ky  (a+1).  ek 


(121) 


(122) 


which  is  equivalent  to  (117). 


COEFFICIENTS  RECURSIVELY  VIA  MOMENTS 


The  starting  point  for  this  case  is  the  characteristic  function  expansion 


in  (101): 

Tf-  m  1_m  ^  m  ^  U+l4m)|.  L. 

f(z)  =  2  C  (-3z)  (1-Bz)  =  ?  c  (-3z)  2  — P - -  (Bz)k  (123) 

m=0  m  m=U  m  k=0  K* 


by  use  of  [5,  15.1.8].  Now  expand  the  left-hand  side  of  (123)  in  powers  of  z, 
according  to  (6),  and  equate  the  coefficients  of  zn  to  get,  for  n  >  0, 


1_ 

nl 


(-*)" 


(a+l+m) 

n-m 

(n-m) l 


n-m 

B 


Cm(-Dm  (an)n 

(n-m) !  (a+l) 

'  '  '  7m 


Therefore 


(o+l)n  Br 


m=0 


nl  cm 
(n-m)i  (a+l)m 


for  n  >  0  , 


(124) 


(125) 
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by  use  of  (92).  Then  using  normalized  moment  definition  (118),  (125)  can  be 
expressed  as 


an  - 


yn  " 


n-1 


m=0 


(-!)'»  (•) 


for  n  >  0 


(126) 


This  is  a  recursive  relation  for  expansion  coefficients  ^an]  in  terms  of 
(normalized)  moments.  The  parameters  a  and  B  in  (118)  are  arbitrary,  except 
that  a  >  -1,  b  >  0. 


SUMMARY 


The  approximations  to  the  probability  density  function  and  cumulative 
distribution  function  are  given  by  (91)  and  (95),  respectively,  where  a  and  b 
are  arbitrary  constants,  except  that  a  >  -1,  b  >  0.  The  generalized  Laguerre 
polynomials  are  available  via  (96).  The  expansion  coefficients  [an]  are 
given  by  the  three  alternatives  (112),  (119),  (126),  in  terms  of  normalized 
cumulants  (62)  and  normalized  moments  (118);  in  the  case  of  (112),  the 
interrelationship  between  expansion  coefficients  [a^  and  [cn]  is  given  in 
(92).  Programs  for  all  three  alternative  procedures  for  determining  the 
expansion  coefficients  [a^  are  listed  in  an  appendix.  The  basis  for  these 
relations  is  the  characteristic  function  expansion  in  (100)-(102). 
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EXAMPLES  OF  HERMITE  EXPANSION 


EXAMPLE  A 


The  first  example  is  one  which  can  be  handled  analytically,  and  thereby 
furnishes  checks  on  numerical  procedures  and  results.  Consider  the  Gaussian 
probability  density  function 

p(u)=-df— )  U  >  0)  ( 127A) 

with  cumulative  distribution  function  and  characteristic  function 

P(u)  -  I  (^)  >  f(1?)  =  exp(ifY  -  \  f2<o2)  .  ( 127B) 

The  cumulants  are 

=  0,  Xi  =  =  u2,  Xn  =  0  for  n  >  3  ,  (128A) 

while  the  moments  are  most  easily  evaluated  by  the  recurrence 

2 

un  =  Y  un_j_  +  ( n— 1)  to  vn_2  for  n  >  2,  =  1,  =  y  •  (128B) 

It  is  obvious  in  this  Hermite  expansion  case  that  the  best  choice  of 
weighting  parameters  would  be  a  =  y>  3  =  to,  for  then  weighting  w  would  match  p 
perfectly  and  there  would  follow  bn  =  0  for  n  2  1*  We  consider  a  mismatched 
choice  of  a  and  0  to  illustrate  rapid  decay  of  the  expansion  coefficients  and 
some  conditions  on  convergence. 

Expansion  coefficient  cn  follows  from  (45)  and  (127A)  according  to 
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cn  =  I'du  p ( u )  =  e  jdx  p(a+ex)  Hen(x)  = 


-s  Hen(x)  He„HS=\,  (129) 


"IVK 


the  last  step  via  use  of  [5,  22.5.18]  and  [8,  7.374  10].  Then  from  (47), 


for  n  >  0  . 


(130) 


This  equation  is  correct  for  all  positive  values  of  8  and  w.  However,  for 
8  <  u,  a  more  convenient  form  can  be  obtained  by  use  of  (74),  if  desired: 


(131) 


,  2  2 
(jo 


where  Hin  is  the  modified  Hermite  polynomial.  For  e  =  u>,  a  limit  of  (130) 


yields  bn  =  (n:)'1/2((Y-a)/6)n. 


If  8  >  u>,  we  can  use  the  result  in  (50A)  on  (130)  and  obtain 

n 


b  or 

n^ 


,2  2 
6  -a) 


3 


-‘A 

n  as  n  >  +  “ 


(132) 


Since  the  quantity  in  parentheses  is  always  less  than  1  in  this  case  of  8  >  u, 
we  have  bn  >  0  as  n  >  +°°. 


For  e  <  id,  we  use  [7,  theorem  8.22.7]  and  find  now  that 


*7  -'/*• 


Yu  -6 


bn  oC  exp(V2n  A)  (  'w  Q~ — )  n  as  n  >  +<*>, 


(133) 
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where  A  is  the  absolute  value  of  the  argument  of  Hen  in  (130).  This 
quantity  (133)  tends  to  zero  with  n,  regardless  of  A,  when  s  >  w/V 21 


Combining  with  the  result  above,  we  can  conclude  that 

bn  >  0  as  n  for  -=l<  b  <  +  oa.  (134) 

n 

Furthermore,  bn  behaves  as  an  n-th  power,  which  is  faster  than  n-^  , 

thereby  guaranteeing  convergence  of  the  probability  density  function  and 
cumulative  distribution  function  series,  according  to  the  discussion  in  (50A) 
et  seq.  On  the  other  hand,  £bn"\  diverges  when  0  <  b  <  w/Y?,  as  may  be  seen 
from  (133). 


The  error  integral  in  (21)  is,  for  Hermite  weighting  (40)  and  probability 
density  function  (127), 


du 


P2(u) 

w(u) 


21 


exp 


,2  2 

2$  -a) 


if 


Y? 


<  B 


(135) 


by  use  of  [8,  3.323  2];  this  integral  is  divergent  if  b  <  w/Y?.  Thus,  for 
this  particular  example,  the  error  integral  and  expansion  coefficient 
sequence  \bn"\  converge  or  diverge  together,  depending  on  the  condition 
B  %  u/Y7.  The  choice  of  a  is  irrelevant  in  this  case. 


A  numerical  example  of  sequence  {b 3  for 

Y  =  1.1,  (U  =  2.3  a  =  1.14,  B  =  2.34  (136) 
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is  plotted  in  figure  1  on  a  logarithmic  ordinate.  Values  of  less  than 
IE— 7  in  absolute  value  are  all  plotted  at  the  ±lE-7  line.  The  critical  ratio 


W /B  in  (130)  is  .184  for  this  example,  leading  to  rapid  decay  of 
expansion  coefficients  {b^ .  The  three  sets  of  expansion  coefficients  in 
figure  1  are  labelled  according  to  the  shorthand  notation 


RC:  Recursively  via  Cumulants, 


DM:  Directly  via  Moments, 


(137) 


RM:  Recursively  via  Moments. 


It  is  seen  that  the  expansion  coefficients  determined  recursively  via 
cumulants,  namely,  the  RC  plot,  decay  rapidly  and  never  encounter  round-off 
error,  whereas  the  DM  and  RM  procedures  both  are  subject  to  large  round-off 
error  for  n  >  70,  as  indicated  by  the  large  increasing  oscillations.  This 
example  can  be  rather  misleading,  however,  since  all  the  cumulants  (128A)  of 
Gaussian  probability  density  function  (127A)  are  zero,  except  for 

o 

=  Y,  =  <o  >  this  leads  to  a  very  special  form  of  the  RC  procedure 
unique  to  the  Gaussian  case. 

In  figure  2,  the  cumulative  distribution  function  and  exceedance 
distribution  function,  l-P(u),  as  determined  by  Hermite  expansion  (48)  using 
N  =  50  terms,  are  plotted.  The  exact  result,  (127B),  overlapped  these  curves 
over  the  full  range  plotted.  The  three  procedures,  RC,  DM,  and  RM,  all 
yielded  identical  distributions  in  figure  2,  as  inspection  of  figure  1 
confirms,  since  the  three  sets  of  expansion  coefficients  are  virtually  the 
same  for  n  <  50.  Even  though  the  three  sets  of  expansion  coefficients  differ 
significantly  for  n  >  60,  the  corresponding  approximate  probability  density 
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n,  sequence  number 

Figure  1.  Hermite  Coefficients  for  Example  fl 


Figure  2.  Distributions  for  Example  fl 
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functions  and  cumulative  distribution  functions  for  N  =  70,  say,  would  not  be 
very  different,  because  the  relative  differences  in  p  and  P  are  very  small, 
somewhere  in  the  IE— 5  range;  see  figure  1  for  n  =  70,  and  recall  that  =  1 
for  this  example. 

EXAMPLE  B 


The  probability  density  function  of  interest  here  is  the  one  previously 
considered  in  (34)  et  seq.: 


p  p 

t  \  2  uY  exp(-u  /a)  )  ^  n 

p(u)  =  y+1  J/La - -  for  u  >  0 


,+1  r(4r) 


(y  >  -1,  id  >  0) 


(138) 


This  class  of  probability  density  functions  includes,  for  y  =  0,  1,  2, 
respectively,  the  one-sided  Gaussian,  Rayleigh,  and  Maxwell  as  special  cases. 
The  characteristic  function  and  cumulants  are  not  easily  determined  directly 
for  this  function.  However,  the  moments,  as  given  already  in  (35),  are 
readily  evaluated  via  the  simple  recursion 


2 

v*n  =  ^-2  *7  ( y- l+n )  Tor  n  >  2,  yQ  -  1, 


(139) 


An  example  of  the  expansion  coefficients  for 


y  =  3,  u=l  a  =  0,  (5  =  .  7 


(140) 


is  depicted  in  figure  3.  The  values  of  bn  for  n  =  0,  1,  2,  3  are  1,  1.90, 
2.18,  1.63,  respectively,  and  lie  above  the  top  of  the  plotted  region.  The 
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i 


Figure  4.  Distributions  for  Example  B 
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coefficients  obtained  directly  via  moments,  DM,  decay  to  approximately  IE— 4 
near  n  =  70  and  then  encounter  round-off  error.  The  expansion  coefficients 
corresponding  to  RC  and  RM  are  more  noisy.  The  procedure  used  for  RC  was  to 
determine  the  moments  via  (139),  transform  directly  to  cumulants  according  to 
(A-7),  and  then  use  (63). 

A  plot  of  the  distributions  using  N  =  65  terms  is  given  in  figure  4;  the 
results  are  the  same  for  all  three  sets  of  expansion  coefficients,  as  may  be 

seen  by  reference  to  figure  3.  Furthermore,  the  exact  cumulative  distribution 

2  2 

function,  P(u)  =  1  -  (l+u  )  exp(-u  )  for  u  >  0,  overlays  these  results 
except  for  the  bow  in  the  exceedance  distribution  function  below  IE— 11  near 
u  =  5.5.  Values  of  the  cumulative  distribution  function  for  u  <  0  as 
determined  by  series  (48)  are  not  zero,  although  they  should  be  for  this 
example;  the  generalized  Laguerre  series  would  fit  this  example  better,  since 
it  is  nonzero  only  for  positive  arguments. 


EXAMPLE  C 


Consider  the  class  of  Bessel-f unction  probability  density  functions 

2  2 

p(u)  =_LuY  exp(-u  lw  )  I^(e-u)  for  u  >  0  , 


(141) 


which  includes  the  Rice  and  generalized  QM  distributions,  for  example.  The 
n-th  moment  is  [8,  6.631  1] 


for  n  >  0  , 


(142) 
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with  h  =  (y+T+l)/2;  in  order  for  mq  to  be  finite,  we  must  have  h  >  0.  The 
^  function  in  (142)  can  be  evaluated  via  recursion;  this  leads  to  a 
recursion  for  the  moments  (see  appendix  E). 


We  consider  here  only  the  special  case  of  the  Rice  probability  density 
function,  namely, 

Jl=  -|  exp  (-  J-J)  ,  y  =  1,  j*=  0  ,  (143) 

0)  ' 

for  which 

P  (u )  =  exp  ~2  -  Ig  (®”U )  for  u  >  0  .  (144) 

0)  \  w  / 

The  moments  in  (142)  then  reduce  to 

un  =  rf?+i)  con  /i  ( §;ii-  j  a2)  ,  d45) 

and  can  be  easily  determined  by  the  recurrence  presented  in  (E-5).  The 
cumulative  distribution  function  corresponding  to  (144)  is  the  Q  function  [1] 


P(u)  =  1 


for  u  >  0  , 


(146) 


the  characteristic  function  is  given  in  [9,  appendix  A]  as  an  infinite  series, 
meaning  that  the  cumulants  cannot  be  determined  directly,  except  via  the 
moments. 


The  particular  example  we  consider  here  for  the  Hermite  expansion  is  a 
sum  of  8  independent  random  variables,  each  with  Rice  probability  density 
function  (144).  For  direct  comparison  with  the  exact  results  in  [9],  we  also 
consider  the  normalized  form  of  (144),  namely  =  2.  Furthermore,  we  limit 
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numerical  consideration  in  this  particular  example  to  evaluation  of  the 
cumulative  and  exceedance  distribution  functions  for  e  =  0,  which  corresponds 
physically  to  the  false  alarm  probability  for  the  sum  of  eight  normalized 
envelopes  of  narrowband  Gaussian  noise  (i.e.,  a  Rayleigh  probability  density 
function  for  the  individual  random  variables). 

For  a  =  4,  e  =  2.15,  the  expansion  coefficients  [b^  are  displayed  in 

figure  5  for  the  RC,  DM,  and  RM  approaches.  All  the  {bn}  for  1  <_  n  <_  20  are 

bigger  than  1;  the  biggest  is  bg  =  12.25.  The  {b  }  ,  for  both  moment 
approaches,  have  not  been  plotted  for  n  >  60  because  they  continue  to 
oscillate  well  beyond  the  ±1  limits,  while  the  RC  coefficients  decay 
exponentially  with  n.  Despite  the  fact  that  the  moments  were  the  initially 
determined  quantities  for  this  example,  the  RC  method  far  outperforms  the  DM 
and  RM  methods,  as  seen  in  figure  5.  The  reason  for  this  is  as  follows:  for 
the  RC  method,  the  procedure  was  to  obtain  moments  via  (145),  cumulants  via 
(A-7),  cumulants  of  the  sum  of  8  independent  random  variables  by  simple 

scaling  by  a  factor  of  8,  and  then  expansion  coefficients  via  (63).  For  the 

DM  and  RM  methods,  the  moments  of  the  sum  of  8  random  variables  were 
determined  via  [6,  (14)]  which  progressively  determined  the  moments  of  a  sum 
of  2  random  variables,  then  3,  4,...,  8  in  order,  and  then  employed  (70). 

This  iterated  procedure  for  moments  requires  more  number-crunching  and  leads 
to  considerably  larger  round-off  error  than  the  simple  scaling  required  for 
the  RC  procedure.  Thus  it  appears  that  when  the  random  variable  of  interest 
is  obtained  as  a  sum  of  several  independent  random  variables,  the  RC  approach 
will  be  the  prime  candidate  for  expansion  coefficient  evaluation;  this  applies 
also  if  the  individual  random  variables  have  different  statistics,  but  remain 
independent. 
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Figure  5.  Hermite  Coefficients  for  Example  C 


Figure  G.  Distributions  for  Example  C 
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The  cumulative  and  exceedance  distribution  functions  for  this  sum  of  8 
normalized  Rayleigh  variates  are  plotted  in  figure  6,  for  the  N  =  140 
expansion  coefficients  of  the  RC  procedure  in  figure  5.  In  order  to  make  a 
precise  determination  of  the  accuracy  of  this  Hermite  series  approach,  the 
false  alarm  probabilities  were  computed  at  the  eight  thresholds  listed  under 
M  =  8  in  [9,  table  1].  To  the  precision  given  in  that  table,  the  computed 
probabilities  were  exactly  the  specified  values  lE-m  for  m  =  1(1)8.  Thus,  as 
anticipated  by  figure  5,  very  accurate  evaluation  of  false  alarm  probabilities 
are  possible  by  this  series  approach. 

A  short  search  of  values  of  the  best  weighting  parameters  a  and  e,  to  use 
with  the  DM  approach,  led  to  a  =  5.84,'  s  =  2.28  and  expansion  coefficients 
bn  near  IE-4  at  n  =  28,  before  round-off  error  became  dominant.  This  is 
better  than  the  result  of  DM  in  figure  5  for  a  =  4,  3  =  2.15.  Evaluation  of 
the  false  alarm  probabilities  at  the  thresholds  in  [9,  table  1]  gave  7  decimal 
accuracy  at  .1,  and  4  decimal  accuracy  at  IE— 8.  This  is  adequate  for  most 
purposes,  but  is  not  as  good  as  the  RC  approach. 

EXAMPLE  D 


In  [4,  appendix  C],  the  characteristic  function  for  shot  noise  with 
random  amplitude  and  duration  modulation,  and  arbitrary  individual  pulse 
shape,  is  derived.  (This  result  is  then  specialized  to  elliptical  pulses  and 
Rayleigh  amplitude  modulation  [4,  (C-36)-(C-42)]. )  Also,  the  cumulants  are 
extracted,  with  general  result  [4,  (24)],  where  v  is  the  average  number  of 
pulses/second,  ^is  the  average  length  of  the  duration  modulation,  ua(n)  is 
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the  n-th  moment  of  the  amplitude  modulation,  and  F(x)  is  the  individual  pulse 
shape  of  the  shot  noise.  Thus  shot  noise  is  a  case  where  the  cumulants  are 
directly  capable  of  evaluation,  whereas  the  moments  must  be  found  indirectly. 

For  the  special  case  of  elliptical  pulses  and  Rayleigh  amplitude 
modulation,  there  follows  for  the  cumulants  [4,  (29)]: 

=  vj  aj  ?  V3(j^l)/f(n+2)  for  n  >  1,  ^  =  0  .  (147) 

These  quantities  are  easily  evaluated  via  recurrence 

=  ^h-2  aa  n2/(n+1)  for  n  >  3,  Xi  =  (jrj 

This  procedure  was  used  in  [4,  appendix  D]  to  obtain  the  probability  density 
function  and  cumulative  distribution  function  results  given  there. 

There  is  a  nuance  that  arises  in  shot  noise  for  pulse  shapes  of  finite 
duration;  see  [4,  pp.  40-42].  Namely,  there  is  an  impulse  in  the  probability 
density  function,  at  u  =  0,  of  area 

p0  =  exp[-vj(x2  -  x1)]  ,  (149) 

where  (x-^,x2)  is  the  non-zero  extent  of  an  unmodulated  individual  pulse. 

Since  an  impulse  is  very  difficult  to  approximate  by  a  finite  series  of 
continuous  functions,  the  effect  of  this  quantity  should  be  subtracted  from 
the  statistics  (moments  or  cumulants),  and  the  continuous  portion  of  the 
probability  density  function  should  be  approximated.  Similarly,  the 
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corresponding  step  in  the  cumulative  distribution  function  at  the  origin 
should  be  eliminated  from  the  approximation  procedure. 

This  feature  is  easily  incorporated  if  Pq  is  subtracted  from  the 
zero-th  order  moment  [4,  p.  42],  The  only  undesireable  side-effect  of  this 
manipulation  is  that  the  initially  computed  cumulants  must  be  transformed  to 
moments,  then  Uq  corrected,  and  then  all  the  new  cumulants  evaluated.  This 
double  transformation  is  necessary  because  the  correction  (subtraction) 
procedure  can  only  be  accomplished  in  the  moment  domain.  Of  course,  when  the 
DM  or  RM  procedures  are  employed  instead  of  RC,  the  last  transformation  to 
cumulants  is  unnecessary;  this  was,  in  fact,  the  procedure  used  in  [4,  p.  60]. 

When  the  individual  pulse  F(x)  has  infinite  duration,  as  for  an 
exponential  or  Gaussian  waveform,  then  X£  -  x^  is  infinite  and  Pq  in 
(149)  is  zero.  In  that  case,  the  considerations  in  the  last  two  paragraphs 
can  be  disregarded,  and  the  cumulants  generated  via  (148)  used  as  is.  It  is 
then  very  likely  that  even  better  accuracy  in  the  expansion  coefficients  will 
be  achieved  than  for  this  current  example. 

For  overlap  factor  [4,  p.  43] 

=  vj(x2  -  x:)  =  6.2,  PQ  =  exp(-6. 2)  =  .00203,  =  1  ,  (150) 

and  for  weighting  parameters  a  =  6.1,  a  =  4.3,  the  expansion  coefficients 
are  displayed  in  figure  7  for  the  three  recursive  procedures.  The  RM 
results  are  considerably  poorer  than  the  RC  and  DM  coefficients,  which  are 
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Figure  7.  Hermite  Coefficients  for  Example  D 
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comparable  for  n  <  65.  However,  even  here,  the  coefficients  have  only  decayed 
to  the  IE— 3  level,  which  may  not  be  sufficiently  small  for  accurate  results. 

The  distributions  using  N  =  65  terms  and  the  RC  expansion  coefficients 
are  given  in  figure  8.  Although  the  actual  cumulative  distribution  function 
is  zero  for  u  <  0,  the  approximation  oscillated  around  zero,  reaching  a 
positive  peak  of  value  .22E-3  at  u  =  -2.  Similarly,  significant  wiggles 
develop  in  the  exceedance  distribution  function  below  the  IE -4  level.  The 
reason  for  the  inadequacy  of  these  Hermite  expansions  near  u  =  0  is  the  abrupt 
zero  behavior  of  the  true  probability  density  function  for  negative  arguments, 
a  feature  inherently  difficult  to  approximate  by  means  of  smooth  continuous 
functions.  The  error  of  the  approximations  in  figure  8  is  estimated  in  a 
later  section  and  superposed  on  the  plot,  for  ease  of  ascertaining  the 
reliability  of  the  curves.  The  corresponding  approximations  for  the 
generalized  Laguerre  series  are  better  for  this  type  of  probability  density 
function,  as  will  be  demonstrated  in  the  next  section. 

The  approximate  probability  density  function  for  this  example,  again  with 
N  =  65  terms,  is  given  in  figure  9  on  a  linear  ordinate.  It  reaches  a 
negative  peak  of  -8E-4,  and  crosses  the  u  =  0  axis  with  value  .004;  both  of 
these  values  should  be  zero,  and  will  be  for  the  generalized  Laguerre  series. 
To  see  how  the  approximate  probability  density  function  behaves  for  larger 
arguments,  the  logarithmic  plot  in  figure  10  is  used.  Wiggles  develop  near 
the  IE-4  level  and  become  large  enough  that  negative  values  of  the  density  are 
yielded  near  u  =  28  and  31.  It  will  be  worthwhile  to  compare  this  Hermite 
series  with  the  generalized  Laguerre  series  to  be  presented  in  the  next 
section.  The  estimated  error  associated  with  figure  10  is  developed  in  a 
later  section. 
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Figure  9.  Linear  Density  for  Example  D 
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EXAMPLES  OF  GENERALIZED  LAGUERRE  EXPANSION 


EXAMPLE  E 


As  with  the  earlier  Hermite  expansions,  the  first  generalized  Laguerre 
example  here  is  one  that  can  be  evaluated  analytically,  for  purposes  of 
checking  numerical  procedures  and  results.  Namely  consider  the  Chi-square 


probability  density  function  of  2(y+1)  degrees  of  freedom  (which  need  not  be 
integer) : 


(151) 


All  probability  density  functions  and  approximations  are  limited  to  u  >  0  in 
this  section,  since  they  are  zero  for  u  <  0;  this  restriction  will  be  presumed 
in  the  remainder  of  the  presentation. 

The  exceedance  distribution  function  is  related  to  the  incomplete  Gamma 
function  [5,  6.5.3]: 


1  -  P(u)  =  dt  p(t)  =  P(y+1,  u/t o)/  T(y+1)  • 


u 


(152) 


The  characteristic  function  follows  from  (151)  as 


f (if)  -  (1  -  ifio)  Y  1 


(153) 


with  cumulants 


(y+1)  n)k  for  k  >  1, 


(154) 


and  moments 


uk  =  (Y+l)k  wk  for  k  _>  0  . 


(155) 
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Thus  either  set  of  statistics  can  be  used  as  a  starting  position.  The  error 
integral  in  (21)  is  finite  if 

-1  <  a  <  2y  +  1  and  e  >  w/2  .  (156) 

We  will  find  the  expansion  coefficients  by  means  of  the  characteristic 
function  expansion  (100),  developed  earlier  for  the  generalized  Laguerre 
series.  Specifically,  we  utilize  the  power  series  expansion 

n  .A-Y-l 


-w/e\ 

.  l-wj  - 

(1-W)Y'“ 

00 

£ 

t— T>. 

m! 

oo 

m  ^ 
w 

m=0 

k=0 

(y+1) 


Wk  , 


(157) 


where  we  used  (153)  and  [5,  15.1.8]  twice.  The  coefficient  of  a  general  term 
wn  is  then  immediately  given  by  the  closed  form 


/s-A"""  „ 

c  =  J>  - ; —  —7 - rr~  f-r— J  for  n  >  0 

n  jjgj  m!  (n-m)l  \  &  / 


Alternative  expressions  for  the  expansion  coefficients  are 
(  y+1), 


(158) 


c  = 
n 


r  i^n  fo-u\n  r  f  B  \ 

sr-  (t7  (a-Y.-m-n-Y;— j  ■ 

(ri)„  /_„An  ,  6n 

-ST- (t)  F 

( a+l)  /  x 

-sr^ F  (-„.n;an;f) 


for  n  >  0  , 


(159) 


obtained  by  means  of  [5,  15.1.1,  15.3.5,  15.3.7]  respectively.  In  fact,  the 
last  result  can  be  obtained  directly  by  using  [8,  7.414  7]  on  (90)  and  (151): 
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oO 


.  _  fH„  uY  exp(-u/a)  ,  (a)  M 

•" '  r  s'  rvu  n  W- 


(160) 


However,  the  latter  two  results  in  (159)  are  not  numerically  stable,  whereas 
(158)  and  the  first  line  of  (159)  are  stable  for  large  n,  without  encountering 
round-off  error. 


Some  special  cases  of  (158)  are  as  follows: 

(y+Dn  /a  \ n 

if  a  =  y,  then  c  = - j —  (- — )  , 

'  n  ni  \  e  / 


(a-t)n 

if  3  =  to,  then  cn  =  — pp —  ; 

if  a  =  y  and  8  =  u,  then  cn  =  5nQ  .  (161) 

The  last  case  is  to  be  expected,  since  the  weighting  exactly  matches  the 
probability  density  function  (151)  then. 


A  numerical  example  of  sequence  (b^  for 

Y  =  1.1,  0)  =  2.3  a  =  1.105,  8  =  2.1  (162) 

is  shown  in  figure  11,  using  the  three  recursive  procedures  developed  earlier 

for  the  generalized  Laguerre  series  in  (112),  (119),  (126).  In  addition, 

exact  result  (158)  is  plotted  for  comparison.  The  expansion  coefficients  have 

a  rapidly  decaying  transient  for  n  <  10,  and  then  a  decay  approximately 
-3  /  2 

proportional  to  n  '  for  large  n.  The  abrupt  change  of  character  at  n  =  5 
does  not  signify  the  onset  of  round-off  error;  rather,  the  latter  is  indicated 
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Figure  11.  Generalized  Laguerre;  Example  E 
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by  an  erratic  behavior,  typically  increasing  exponentially  with  n  (linear 
growth  on  a  logarithmic  ordinate). 

A  different  plotting  strategy  will  be  adopted  henceforth  for  the 
expansion  coefficients  {bnj ,  in  order  not  to  clutter  the  diagrams  with  large 
oscillations  as  in  figures  1,  3,  5,  7.  Specifically,  when  the  expansion 
coefficient  bn  first  exceeds  the  ±1  limits,  the  remainder  of  sequence  [bn] 
will  not  be  plotted,  since  this  is  a  region  of  large  round-off  error.  Thus, 
although  the  RM  curve  in  figure  11  returns  to  the  =*=1  limits  briefly  at 
n  =  52,53,  these  values  are  not  displayed. 

Round-off  error  for  the  RC  procedure  does  not  become  as  significant  as 
for  the  two  moment  approaches  until  n  has  increased  by  almost  10,  for  this 
example  in  figure  11.  In  fact,  the  expansion  coefficients  for  the  RC 
procedure  overlap  the  exact  values  until  n  =  40.  The  corresponding 
approximate  distributions,  using  N  =  40  terms  in  expansion  (95)  as  determined 
by  RC,  are  plotted  in  figure  12.  The  exact  result  (152)  overlays  these 
results  over  the  entire  range  plotted. 


EXAMPLE  F 


The  following  probability  density  function  corresponds  to  a  noncentral 
Chi-square  variate  of  2v  degrees  of  freedom: 


Iv_  i(dVT?) 


(v  >  0)  ; 


(163) 
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d  is  the  noncentrality  parameter,  and  2v  need  not  be  integer, 
characteristic  function  is  [8,  6.631  4] 


The 


f(i^)  =  ( 1-i 2f )  'v  exp(^j|j)  . 


(164) 


and  is  the  same  as  the  one  considered  in  [10,  (50)  et  seq.].  The  exceedance 
distribution  function  is  the  generalized  Q-function: 

1  -  P(u)  =  dt  \  exp  (-  (dYt1)  = 


=  dx  x  exp  (-  -  ^  j  (j)  Iv-1(dx)  =  Qv(d,YuO  • 

By  expanding  the  /Pn  of  (164)  in  a  power  series  in  if,  the  cumulants 
follow  as 

=  2n( n-1).'  +  j  d2  n^  for  n  >  1,  ^  =  0  . 


(165) 


(166) 


And  the  moments  are  obtained  from  (163)  as 


v»n  =  2n  (v)n  ^F^(-n; v,-d?/2)  = 

=  2n  n!  (-d2/2)  for  n  >  0  ,  (167) 

by  use  of  [8,  6.631  1]  and  [5,  13.6.9].  Both  (166)  and  (167)  lend  themselves 
to  simple  recurrences  which  involve  only  positive  quantities;  thus  the 
starting  statistics  can  be  quickly  and  accurately  evaluated. 


The  numerical  example  we  consider  here  will  be  compared  with  the  exact 
results  in  [10,  figure  11],  namely, 
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v  =  2.7,  d  =  3  a  =  1.7,  £  =  5.5  .  (168) 

Since  the  probability  density  function  in  (163)  behaves  as  uv--*-  as  u  >  0+, 
it  is  reasonable  to  choose  weighting  parameter  a  in  (82)  as  v-1,  as  indicated 
in  (168).  And  since  (163)  behaves  as  exp(-u/2)  as  u  +o6,  we  must  choose 
S  >  1  in  order  that  the  error  integral  in  (21)  is  finite.  The  particular 
values  in  (168)  approximately  minimize  the  sum  of  {b^jg  in  (21). 

The  expansion  coefficients  {bn"[  as  determined  by  the  three  available 
recursive  procedures  are  displayed  in  figure  13.  The  RC  coefficients  decrease 
to  values  less  than  IE— 10  near  n  =  50,  before  round-off  error  becomes 
significant.  The  two  moment  approaches  deteriorate  near  n  =  30,  which  is 
markedly  poorer  than  the  cumulant  approach.  The  distributions,  as  determined 
by  N  =  50  terms  of  the  RC  approach,  are  given  in  figure  14,  and  agree  with  the 
d  =  3  curve  of  [10,  figure  11].  When  the  approximate  probability  density 
function  for  N  =  50  was  compared  with  exact  result  (163),  10  decimals  of 
agreement  were  obtained;  this  is  due  to  the  ability  to  get  very  small  ^b  | 
in  figure  13  via  the  RC  method. 

EXAMPLE  G 


This  example  is  the  Rice  probability  density  function  given  in  (144), 
with  moments  (145)  and  cumulative  distribution  function  (146).  The  starting 
statistics  are  the  moments  as  determined  by  recurrence  ( E— 5) — ( E— 6) . 
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Figure  13.  Generalized  Laguerre;  Example  F 
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Figure  14.  Distributions  for  Example  F 
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The  particular  numerical  case  of  interest  is 

e  =  3,  o)2  =  2  a  =  1,  6  =  1.  (169) 

The  values  of  a  and  e  were  found  by  the  usual  trial  and  error  search  procedure 
of  observing  plots  of  expansion  coefficients  {bn} ,  looking  for  rapid  decay 
and  small  round-off  error;  results  for  this  example  are  displayed  in  figure 
15.  The  RM  procedure  deteriorates  rapidly  at  n  =  30,  whereas  DM  and  RC  are 
useable  up  to  n  =  55  and  65  approximately. 

The  cumulative  and  exceedance  distribution  functions  for  N  =  65  terms  of 
the  RC  procedure  are  plotted  in  figure  16,  along  with  exact  result  (146).  The 
approximate  exceedance  distribution  function  overlaps  the  exact  one  until 
slightly  below  the  probability  level  IE-4,  which  corresponds  to  the  level  of 
reliability  of  bn  in  figure  15  at  n  =  65.  Then  the  exceedance  distribution 
function  makes  a  positive  (upward)  turn  below  IE -6,  which  is  impossible  for  a 
physical  density  function  which  must  remain  positive;  thus  the  approximation 
deteriorates  rapidly  for  u  >  7. 

EXAMPLE  H 

This  is  a  follow-on  to  the  previous  example,  in  that  we  consider  a  sum  of 
8  Rice  variates,  each  with  the  statistics  in  (169).  The  expansion 
coefficients  for 

e  =  3,  o)2  =  2  a  =  26,  6  =  1  (  170) 
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n,  sequence  number 

Figure  15.  Generalized  Laguerre;  Example  G 


Figure  16.  Distributions  for  Example  G 
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are  displayed  in  figure  17.  Whereas  both  DM  and  RM  are  useless  beyond  n  =  25, 
the  expansion  coefficients  determined  by  RC  decay  down  to  the  IE— 10  level  at 
n  sr  150  before  round-off  error  becomes  significant.  The  corresponding 
distribution's  in  figure  18,  using  N  =  143  terms  of  the  expansion  via  RC, 
reveal  accurate  results  down  to  the  IE— 12  level  of  probability,  except  for  a 
slight  flare  in  the  exceedance  distribution  function  below  IE— 11. 

We  also  checked  the  example  of  the  sum  of  8  normalized  Rayleigh  variates 
considered  earlier  via  a  Hermite  series  in  example  C.  For  a  =10,  $  =  .9,  the 
expansion  coefficients  {bn}  decayed  to  the  IE— 11  level  at  n  =  100  for  the  RC 
approach  and  agreed  with  the  false  alarm  probabilities  calculated  exactly  in 
[9,  table  1]  for  M  =  8.  By  contrast,  the  DM  expansion  coefficients  were 
subject  to  significant  round-off  error  by  the  time  n  reached  30,  and  were 
useless  for  small  probability  calculations. 

EXAMPLE  I 

We  return  to  the  shot  noise  process  previously  considered  via  a  Hermite 
series  in  example  D.  The  equations  and  discussions  there  should  be  reviewed, 
since  they  are  directly  relevant  to  the  generalized  Laguerre  expansion  here. 

For  the  choice  of  parameters  in  (150),  the  selection  of  generalized  Laguerre 
weighting  parameters 

a  =  .74,  6  =  2. 1  (171) 

leads  to  the  expansion  coefficients  plotted  in  figure  19.  The  DM  and  RC 
results  agree  to  n  =  32,  and  then  begin  to  diverge  from  each  other.  By  way  of 
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Figure  17.  Generalized  Laguerre;  Example  H 


Figure 


u ,  th r esho  1  d 

18.  Distributions  for  Example  H 
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n,  sequence  number 

Figure  19.  Generalized  Laguerre;  Example  I 


Figure  20.  Distributions  for  Example  I 
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contrast  with  the  Hermite  expansion  coefficients  in  figure  7,  where  values  in 
the  IE— 3  range  were  achieved,  values  in  the  IE— 6  range  can  be  obtained  here 
for  the  generalized  Laguerre  expansion,  for  n  in  the  rni d— 30s -  The  DM  result 
was  previously  given  in  [4,  figure  D-l]. 

The  distributions  for  M  =32  terms  of  the  RC  procedure  are  plotted  in 
figure  20.  This  result  is  considerably  better  than  the  Hermite  expansion  in 
figure  8;  instead  of  the  wiggles  which  developed  at  IE— 4  in  figure  8,  the 
curve  in  figure  20  is  smooth  down  to  the  IE -8  probability  level,  and  then 
develops  a  bump.  Also,  the  cumulative  distribution  function  is  accurate  at 
u  =  0,  where  it  takes  on  the  value  PQ  =  .002  given  in  (150),  and  is  zero  for 
u  <  0.  This  cumulative  distribution  function  was  previously  given  in 
[4,  figure  3]. 

The  probability  density  function  for  N  =  32  terms  of  the  RC  procedure  is 
given  in  figure  21;  this  result  was  previously  given  in  [4,  figure  9].  It  is 
significantly  better  near  the  origin  than  the  Hermite  approximation  given 
earlier  in  figure  9,  which  developed  negative  values  for  u  <  0.  In  order  to 
see  what  the  probability  density  function  does  for  larger  u  values,  the  same 
probability  density  function  is  plotted  on  a  logarithmic  ordinate  in  figure 
22.  It  is  accurate  to  the  IE  —9  level  but  then  develops  a  hook  that  is 
incorrect;  however,  this  approximation  remains  positive  even  at  this  very  low 
value  of  the  density,  whereas  the  corresponding  result  via  a  Hermite  expansion 
in  figure  10  developed  negative  values.  The  estimated  errors  in  figures  20 
and  22  are  evaluated  in  a  later  section. 
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EXAMPLE  J 


This  last  example  is  for  probability  density  function 
p  (u )  =  ^  exp  (h/1)  for  u  >  0  , 


(172) 


for  which  the  moments  are 


un  =  (2n+l)l  . 


(173) 


The  characteristic  function  and  cumulants  are  not  available  in  any  convenient 
analytic  form. 

This  is  a  particularly  difficult  example,  since  the  characteristic 
function  expansion  in  (6)  has  a  zero  radius  of  convergence;  thus  the  moments 
do  not  uniquely  determine  the  probability  density  function  or  cumulative 
distribution  function.  Also,  the  error  integral  in  (21)  is  always  infinite; 
in  fact,  regardless  of  the  choice  of  weighting  parameters  a  and  @  used  in  the 
generalized  Laguerre  series,  the  expansion  coefficients  ^b  ^  always 
diverged.  Nevertheless,  a  search  of  parameter  values  led  to  a  pair  of 
selections,  namely, 


(174) 


a  =  -. 35,  0  =  30, 


for  which  the  expansion  coefficients  had  an  initial  decay  to  the  IE— 2  level 
before  divergence  took  over;  see  figure  23.  In  fact,  the  identical  same 
results  were  obtained  for  all  three  methods,  RC,  DM,  RM;  this  is  probably  due 


to  the  fact  that  divergence 
significant. 


dominated  before  round-off  error  became 
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n,  sequence  number 

Figure  23.  Generalized  Laguerre;  Example  J 


Figure  24.  Distributions  for  Ex  amp  1 e  J 
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The  distributions  are  plotted  in  figure  24  for  N  =  15  terms  of  the 
generalized  Laguerre  series.  Comparison  with  the  exact  exceedance 
distribution  function 

1  -  P(u)  =  (1  +  u1/2)  exp  (-ui/2)  for  u  >  0  (175) 

reveals  that  the  approximation  is  decent  down  to  the  .01  probability  level, 
but  then  oscillates  more  and  more  violently  as  u  increases.  Thus  even  in  this 
non-unique  example,  a  limited-quality  approximation  is  achieved  by  the 
generalized  Laguerre  series;  this  example  confirms  the  comment  in  [3,  p.  167] 
that,  even  for  a  divergent  series,  a  limited  number  of  expansion  coefficients 
often  gives  a  satisfactory  approximation. 

The  exact  and  approximate  probability  density  functions  are  plotted  on  a 
linear  ordinate  in  figure  25,  and  on  a  logarithmic  ordinate  in  figure  26, 
using  N  =  15  terms  of  the  generalized  Laguerre  series,  when  the  expansion 
coefficients  were  determined  by  the  DM  method.  The  approximate  probability 
density  function  is  negative  for  150  <  u  <  190,  around  the  IE —6  level.  The 
estimated  errors  of  the  approximations  in  figures  24  and  26  will  be  developed 
in  the  next  section. 
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Figure  25. 


Linear  Density  for  Example  J 


Figure  26.  Log  Density  for  Example  J 
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ESTIMATED  ERRORS  OF  APPROXIMATIONS 


When  the  calculations  of  the  approximate  cumulative  or  exceedance 
distribution  functions  or  the  corresponding  probability  density  function  are 
made,  it  would  be  very  useful  to  have  a  rough  estimate  of  their  reliability. 
One  way,  as  discussed  in  the  previous  sections,  is  to  look  for  nonsmooth  or 
anomalous  behavior  on  the  tails  of  the  functions.  Here,  we  will  develop  a 
more  quantitative  estimate  of  the  error  and  superpose  it  on  some  of  the 
previous  examples,  for  confirmation. 

Both  the  Hermite  and  generalized  Laguerre  orthonormal  polynomials 
oscillate  with  n  and  decay  slowly.  The  same  general  behavior  is  true  of 
expansion  coefficients  [bn|.  This  leads  to  summations  for  the  various 
functions  with  terms  that  also  oscillate  and  decay.  A  rough  estimate  of  the 
error  is  afforded  by  the  envelope  of  these  oscillations,  evaluated  at  the 
first  neglected  term  of  the  summation.  This  procedure  will  be  pursued  for 
both  types  of  expansions;  how  useful  it  is  will  be  indicated  by  numerical 
examples. 

HERMITE  EXPANSION 

The  following  result  for  the  envelope  of  the  Hermite  polynomial  is 
obtained  from  [5,  6.1.39  and  22.5.18]  and  [7,  8.22.8]: 


(176) 


Also,  from  (46)  and  (47),  the  n-th  term  of  the  approximate  probability  density 
function  is 
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i b^K  Hs(r -)  • 


(177) 


Then  the  magnitude  of  the  error  of  the  probability  density  function 
approximation,  if  the  n-th  term  is  the  first  one  neglected,  is  roughly 


E„  a 


Envjbj  Env[(ni)  A  Hen(^)]  - 


?'A  * 

L  ir  Q 
L 


exp  /-  a-l;  -J  n  Env[b^  as  n  >  +c<5, 


2\  -»/« 


(178) 


Here  we  used  (176) . 


As  for  the  cumulative  distribution  function,  we  have  from  (47)-(49),  the 
n-th  term  of  the  approximation  as 

-<»(¥)  bn  Hen-1  (*?)■  <179> 


The  magnitude  of  the  error  for  the  cumulative  and  exceedance  distribution 
functions,  if  the  n-th  term  is  the  first  one  neglected,  is  then  defined  as 


En(U;p)  »  i  Env  {bj  Env  [(n!)  A  He^  (■“=*)]  - 

~  [2 ^ir  exp  ^  y  n  Env  [bn"j  as  n  >  +M.  ( 180) 


Again,  (176)  was  of  crucial  importance  in  getting  this  result. 
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Since  the  above  estimates  are  asymptotic  in  n,  they  will  be  most  reliable 
for  n  large;  their  use  for  small  n  could  be  very  misleading.  The  way  to  use 
these  error  estimates  for  the  density  and  distribution  approximations  is  as 
follows.  First,  a  search  on  a  and  e,  to  find  the  fastest  decaying  expansion 
coefficients  [bj ,  is  conducted.  The  weighting  parameter  values,  a  and  e, 
and  the  corresponding  envelope  value  of  the  expansion  coefficients  {bn}  at 
the  point,  n,  where  round-off  error  becomes  dominant,  are  then  noted.  (For 
example,  for  figure  7,  we  observe  that  Envlb^  z  2E-3  at  n  =  65,  when 
a  =  6.1,  3  =  4.3;  see  example  D.)  Then  (178)  and  (180)  can  be  computed  and 
plotted  in  the  ranges  of  u  of  interest. 

An  example  of  this  procedure  for  the  shot  noise  process  in  example  D  is 
given  in  figures  27  and  28.  In  particular,  the  approximate  results  are 
repeated  from  figures  8  and  10,  and  error  measures  (180)  and  (178), 
respectively,  are  superposed  as  dashed  lines,  each  on  the  appropriate  figure. 
Just  where  the  approximations  develop  large  wiggles,  the  errors  are  of 
comparable  magnitude,  indicating  unreliable  estimates  there. 

It  should  be  observed  from  these  figures  (or  from  (178)  and  (180))  that 
the  absolute  error  is  maximum  at  u  =  a,  but  that  the  relative  error  is  a 
minimum  in  that  neighborhood.  Also,  although  the  absolute  error  decays  with 
u,  the  correct  answer  decays  faster,  leading  to  an  increasing  relative  error, 
which  eventually  becomes  so  excessive  in  the  tails  of  the  various  functions 
that  the  approximations  are  useless. 
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Figure  28.  Estimated  Error  of  Figure  10 
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GENERALIZED  LAGUERRE  EXPANSION 


The  details  for  the  generalized  Laguerre  series  are  very  similar  to  those 
above  and  so  will  be  abbreviated.  The  envelope  of  the  generalized  Laguerre 
polynomial  is  [7,  8.22.1] 

*  _1  X_  _a  _1  _a  _1 

Env^L^“^(x)|  -  ir^e^x^  ^  n^  ^  as  n  >  +  ot>,  for  x  >  0  .  (181) 


From  (91)  and  (92),  the  n-th  term  of  the  approximate  probability  density 
function  is 


ua  exp(-u/s) 

ea+1  ru+l) 


(a) 


(182) 


Then  the  magnitude  of  the  error  of  the  probability  density  function 
approximation  is,  for  u  >  0, 

E  (u;p)  =  Env{b1  Env((Ufj-Y*L<»)  - 

n  6“+i  p(<»+i)  nJ  w°nv  " 

T  1'/2 

- Wrto.+l)e2j  (jj  exp(-  2b)  n  Env  Ibnlas  n  »  +co> 


(183) 


where  we  used  [5,  6.1.47]  and  (181).  This  quantity  peaks  at  u  =  0(01-5-). 


With  regards  to  the  cumulative  distribution  function,  the  n-th  term  of 
the  approximation  is,  from  (95)  and  (92), 
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+  1 

u  exp (-u/e) 

e“+1  rc»-D 


(184) 


Then  the  magnitude  of  the  distribution  error,  for  both  the  cumulative  and  the 
exceedance  distribution  functions,  is  roughly 


as  n  >  +  *>,  for  u  >  0,  (185) 


upon  use  of  (181).  This  quantity  reaches  its  peak  at  u  =  s(a+i). 

An  application  of  these  results  to  the  shot  noise  process,  which  was 
re-investigated  in  example  I  via  the  generalized  Laguerre  series,  is  given  in 
figures  29  and  30.  Specifically,  the  approximate  results  from  figures  20  and 
22  have  been  repeated,  and  error  measures  (185)  and  (183),  respectively, 
superposed  as  dashed  lines.  They  confirm  the  earlier  observations  that  the 
distribution  and  density  approximations  are  reliable  until  the  anomalous 
behavior  on  the  tails  manifests  itself. 

The  difficult  example  J  is  considered  in  figures  31  and  32.  Since  the 
expansion  coefficient  sequence  {bn^  in  figure  23  diverged  for  large  n,  the 
selection  of  n  =  15,  as  used  in  figures  24-26,  is  not  the  large  value  needed 
to  justify  the  use  of  (183)  and  (185).  Thus,  the  dashed  curves  on  figures  31 
and  32  must  be  considered  only  as  ball-park  estimates;  in  general,  the 
approximate  error  appears  to  be  too  conservative  in  these  two  figures. 
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Figure  29.  Estimated  Error  of  Figure  20 


u,  threshold 

Figure  30.  Estimated  Error  of  Figure  22 
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Figure  31.  Estimated  Error  o-f  Figure  24 


Figure  32.  Estimated  Error  of  Figure  26 
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Finally,  the  Rice  variate  of  example  G  is  re-considered  in  figure  33.  We 


extrapolating  in  figure  15  from  smaller  n. 


since  round-off  error  is  becoming  significant  by  this  value.  It  verifies  the 
unreliability  of  the  approximation  in  figure  33  for  u  >  7. 

Although  all  the  examples  in  this  report  have  the  capability  of 
evaluating  either  the  moments  or  the  cumulants  via  recursion,  this  is  by  no 
means  necessary.  Any  method  whatsoever  of  accurately  calculating  the  starting 
statistics,  be  they  moments  or  cumulants,  is  acceptable.  For  example,  if  a 
random  variable  with  known  probability  density  function  q  is  passed  through  a 
complicated  nonlinearity  g,  the  moments  of  the  output  are  given  by 


(186) 


These  quantities  could  be  evaluated  for  0  <_  n  <_  N  by  brute-force  numerical 
procedures  if  necessary.  The  limit  value  N  will  depend  on  the  accuracy  with 
which  g  and  q  can  be  evaluated,  if  g(u)  _>  0  for  all  u,  these  integrals  can  be 
accomplished  to  a  high  degree  of  accuracy,  thereby  allowing  large  values  of  N 
to  be  employed. 
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Figure  33.  Estimated  Error  of  Figure  16 
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DISCUSSION 


Several  alternative  methods  have  been  presented  for  obtaining  either 
Hermite  or  generalized  Laguerre  series  expansions  of  probability  density 
functions  or  cumulative  and  exceedance  distribution  functions,  by  means  of 
recursive  relations  involving  either  moments  or  cumulants.  Furthermore, 
estimates  of  the  errors  of  the  approximations  are  furnished  so  that  the 
reliability  can  be  assessed.  Comparisons  between  approximations  obtained  by 
either  the  Hermite  or  the  generalized  Laguerre  series  afford  an  assessment  of 
the  accuracy  of  each;  also,  the  availability  of  three  alternative  recursive 
procedures  for  the  expansion  coefficients  allows  for  selection  of  the  best 
method  and  results,  and  determination  of  the  amount  of  round-off  error. 

The  key  feature  to  this  approach  is  the  rapid  calculation  and  observation 
of  the  orthonormal  expansion  coefficients  |b^  for  each  particular  guess  of 
weighting  parameters  a  and  e.  A  trial  and  error  procedure  is  suggested  for 
determining  o  and  s  values  that  yield  the  set  of  fastest-decaying  expansion 
coefficients.  From  observation  of  the  expansion  coefficients,  the  number  of 
terms  to  retain  in  the  series  expansions  is  ascertained,  being  sure  to  avoid 
the  effects  of  round-off  error  which  dominates  the  calculated  expansion 
coefficients  (b^  for  large  n.  Since  the  amount  and  location  of  round-off 
error  on  the  plot  of  expansion  coefficients  also  depends  on  a  and  e,  a 
judicious  search  may  be  required  to  find  acceptable  weighting  parameter 
values.  Of  course,  a  computer  with  a  larger  number  of  significant  digits 
would  greatly  alleviate  this  drawback;  the  particular  computer  used  for  all 
the  calculations  reported  here  is  the  Hewlett-Packard  9000  Model  520  which 
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devotes  52  bits  (15.65  decimal  digits)  to  the  mantissa  and  11  bits  to  the 
exponent.  Failure  of  the  technique  is  indicated  by  divergence  of  the 
expansion  coefficient  sequence  [b^. 

Programs  for  the  shot  noise  process  considered  in_examples  D  and  I  are 
presented  in  appendix  F.  Times  of  execution  are  as  follows.  For  the  Hermite 
series,  the  80  cumulants  or  80  moments  required  as  input  for  figure  7  took  .7 
or  .35  seconds,  respectively.  The  calculation,  plotting,  and  display  of  the 
80  expansion  coefficients  in  figure  7  took  1.6  seconds  via  the  RC  approach  and 
1.75  seconds  via  the  two  moment  approaches.  The  computation  and  display  of 
the  100-point  plots  of  the  cumulative  distribution  function  in  figure  8  and 
the  probability  density  function  in  figure  9,  each  using  65  terms  in  the 
series  expansion,  took  1.1  and  .95  seconds,  respectively. 

For  the  generalized  Laguerre  series,  the  70  cumulants  or  70  moments 
required  as  input  for  figure  19  took  .54  seconds  or  .28  seconds, 
respectively.  The  calculation  and  display  of  the  70  expansion  coefficients  in 
figure  19  took  1.8  seconds  via  the  RC  approach  and  1.5  seconds  via  the  two 
moment  approaches.  The  computation  and  display  of  the  100-point  plots  of  the 
cumulative  distribution  function  in  figure  20  and  the  probability  density 
function  in  figure  21  took  1.1  and  .7  seconds,  respectively.  These  execution 
times  are  short  enough  to  allow  a  human  observer  to  conduct  a  rapid 
trial -and-error  search  of  a,e  space,  determine  adequate  parameter  values,  and 
assess  their  accuracy. 

Alternative  exact  procedures  for  determination  of  cumulative  and 
exceedance  distribution  functions  via  characteristic  functions  have  been 
presented  in  [9,  10,  11].  Those  methods  generally  have  the  potential  for 
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greater  accuracy,  are  less  subject  to  round-off-error,  and  would  be  preferred 
if  possible.  However,  analysis  of  systems  with  nonlinearities  and  memory 
sometimes  precludes  or  greatly  hinders  their  application;  in  such  cases,  the 
current  approach  is  a  very  good  candidate  for  consideration. 

The  two  weightings  in  (1)  and  (2),  namely  the  Hermite  and  generalized 
Laguerre,  have  been  investigated  rather  intensively  here,  because  so  many 
properties  and  recursions  are  available  for  the  corresponding  (orthonormal) 
polynomials.  These  properties  have  been  utilized  to  derive  simple  recursive 
relations  for  the  expansion  coefficients  and  density  and  distribution 
functions,  thereby  realizing  quick  efficient  procedures  for  numerical 
evaluation  and  observation. 

It  would  be  extremely  useful  to  be  able  to  extend  these  results  to  the 
weighting 

9  9 

u“  exp(-u  /e  )  for  u  >  0,  (187) 

since  this  class  of  probability  density  functions  is  often  encountered  in 
nonlinear  systems  with  Gaussian  inputs.  However,  there  are  several  pivotal 
recursive  relations  for  the  corresponding  orthonormal  polynomials  that  would 
be  needed,  and  it  is  questionable  if  a  fast  procedure  could  be  devised  without 
them.  Also,  it  is  unknown  if  recursive  procedures  for  the  expansion 
coefficients  in  terms  of  moments  or  cumulants  could  be  derived,  as  was  done 
here  for  the  Hermite  and  generalized  Laguerre  weightings.  This  is  a  topic 
worthy  of  further  investigation. 
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APPENDIX  A.  COEFFICIENT  RECURSION  FOR  EXPONENTIAL  OF  POWER  SERIES 


oo 


Suppose  power  series  hn  z11  converges  for  some  lz|  >  0,  and  we 

n=0  n 


exponentiate  it,  getting  a  new  power  series 
oo  r  oo  I 

2  9nzn»exp|Z 

n=0  k  n=0  j 


Then  the  lowest  order  coefficient  is 

g0  =  exp(hQ)  , 

while  for  k  _>  1,  we  have 

%  =YT  (af) 


( A— 1) 


(A— 2) 


i 

r  oo 

exp  1 

2  h„z4 

Ln=0  J 

-*z=0 


1  f  d^f  f  h  n-l  ( 5  . 

ln=inhnZ  8XPlShnZ 


ki  l^dz 

*  %  m 


n 


n  h  z 


n-l 


n=l 


'z=0 


d 

dz 


z=0 

k-l-p 


o c 


exp"S^_ .  h  z 
(n=0  n 


z=0 


-IT  2  (V)  (p+1);  hpH  (k-!-P)-  9k_l_p  - 


k-1 

■c 

P 
k-1 


-  k  ^  ^  ^p+l  gk-l-p  “  k  gk-m  * 

p=0  r  r  m=l 

Thus  we  have  the  recursion  for  coefficients  fg^  in  terms  of  the  £hm^: 

1  J- 

gk  =  k  ^  m  hm  gk-m  for  k  1  g0  exP(h0}  * 

m=l 


( A— 3 ) 


( A— 4) 
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If  we  now  refer  to  (6)  and  (7)  and  identify 

g„-Pn/n:-  hn=Vn:- 

there  follows  the  moments  in  terms  of  the  cumulants  according  to 


(A-5) 


uk  = 


for  k  >  1, 


u0  =  expQy  . 


( A  —6 ) 


This  is  a  slight  generalization  of  [6,  (10)].  This  equation  is  immediately 
inverted,  to  yield  the  cumulants  in  terms  of  moments: 


X  = 
k  u0 


for  k  >  1, 


-An  Pq  , 


( A  —7 ) 


which  generalizes  [6,  (11)]. 


In  terms  of  the  normalized  cumulants  and  moments  defined  in  (62)  and 
(69)  respectively,  we  have 
k-1  A 


A  1 

uk  =  k 


m=0 


•^k-m  ym  for  k  ^  *0  =  exP(^0} 


( A-8) 


(A -9) 
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APPENDIX  B.  EXPANSION  OF  Hen(x+y) 


The  quantity  Hen(x+y)  is  a  polynomial  of  degree  n  in  y.  Therefore  we 


can  expand 


_ _  „m 

He  (x+y)  =  S  y  i  > 
n  m  ml  ’ 


( B— 1 ) 


where  Ym  will  also  depend  on  n  and  x.  In  fact. 


^  =  (l7)ra[Hen(X+,)]y,0"(i)m[H%(t»]t=x- 

"(it)  [n  HVl(t)lt_x  -  n(n-1>  •••  H,WX>  * 

-  iSisTT  Hen-m^x^  • 

where  we  used  [5,  22.8.8]  repeatedly.  Using  (B-2)  in  (B-l),  we  have  the 
alternative  forms  for  the  expansion. 


(B-2) 


Hen(x+y)  = 


|  (:)  ?  - 

5  fn 

,=o 

|  (c)  x 


x  - 


for  n  >  0 


(B-3) 
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APPENDIX  C.  EVALUATION  OF  In(y)  IN  (94) 


We  have,  from  (94), 


I  (y)  = 

nVJ  / 


dx  xa  e  x  L^(x)  for  n  ^  0  . 


(C-l) 


Then 


pJ  a+l  -y 

I0(y)  =  dx  xa  e  x  1  =  y(a+l,y)  =  y-  a+|-e-  ■  XF ^Ua+Ziy) t  (C-2) 

using  [5,  22.4.7,  6.5.2,  and  6.5.12].  Also,  we  have  from  [5,  22.11.6], 

x“  e-x  !>>(x)  -  4  Ml"  ?e‘x  x'**"]  .  (C-3) 

n  nl  \dx/  i  J  v  ‘ 

Then  for  n  ^  1,  (C-l)  can  be  developed  as 


,n-l 


>-y  y“+n]  - 


(C-4) 


where  we  set  the  lower  limit  of  the  evaluated  integral  to  zero  since 
a+n  >  a+l  >  0. 
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APPENDIX  D.  FOURIER  TRANSFORM  OF  GENERALIZED  LAGUERRE  POLYNOMIAL 


We  wish  to  evaluate  transform 
oo 

A(.)  -  J  dt  e1"*  t“  e-t  L(“>(t)  . 


Now 


n!  t“  e-‘  L<*>(t)  =  (jf)"  [e-‘  for  n  >  0 

according  to  [5,  22.11.6].  Therefore  for  n  >  1, 

n:  AM  =.  J  dt  e'“l  £e-t  t“  "}  ■ 

-  ^ ta+n]  - 
■  r « ^  (d Cl  {•-*  ‘“i . 


(D-l) 


( D— 2 ) 


( D— 3) 


where  we  used  integration  by  parts  with  the  fact  that  the  integrated  part  is 
zero  at  t  =  0  and  >o  ,  since  a+n  a+l  >  0.  Repeated  integration  by  parts  then 
yields 

n!  A( oj)  =  (-iu)n  [  dt  eiut  e_t  ta+n  =  p(a+l+n)  .  (D-4) 

o  (l-iu)a 

This  is  the  result  quoted  in  (104). 
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APPENDIX  E.  RECURRENCE  FOR  EXAMPLE  C 


The  starting  point  is  the  moment  expression  in  (142): 
*  p(|+h)  .  n+2h 


*  el  \T"  «"  r  /niL  ^  \ 

Un  -  lFirh’W’Z)  • 

where  h  =  (y+J+l)/2,  z  =  u>2e2/4.  Denote  the  ^  term  in  (E-l)  by 
F  ,  and  the 
immediately 


Fn,  and  the  leading  factor  by  Gn;  thus  un  =  Gn  F  .  There  follows 


G  =  G  i»‘ 
n  n-2 


!(M 


for  n  >  2 


For  the  ]F]_  function,  we  refer  to  [5,  13.4.1]  to  get 
Fn  ■  a ^7  [ln*2h-3-rz)Fn.2  +  (*♦ 2-  £  -l)  Fn_4]  . 


(E-l) 


(E-2) 


( E-3) 


If  we  substitute  (E-2)  and  (E-3)  into  un  =  Gn  Fn,  and  then  re-apply 
(E-2)  in  the  second  term,  we  obtain 

un  =  u)2(n+Y-2+z)  un_2  -  ^  a)4  ^(n+y-3)2  -[f2]  un_4  ;  (E-4) 

we  also  eliminated  h.  Starting  values  for  yn  can  be  obtained  from  (E-l). 


For  the  special  case  (143)  and  (144),  (E-4)  reduces  to 

2,  v  1  4,  2 

un  =  to  (n-l+z)  un_2  -  4  »  (n-2)  un_4  , 


(E-5) 
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with  starting  values 


(E-6) 


Kummer's  transformation  [5,  13.1.27]  was  employed  in  this  last  equation;  these 
forms  afford  accurate  starting  values  for  recursion  (E-5). 
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APPENDIX  F.  PROGRAM  LISTINGS 

Eight  programs  are  listed  in  this  appendix.  They  are  given  in  BASIC  for 
the  Hewlett  Packard  9000  Model  520  computer.  For  ease  of  reference,  a 
shorthand  notation  is  adopted: 


p 

denotes 

cumulative  or  exceedance  distribution  function 

p 

denotes 

probability  density  function 

H 

denotes 

Hermite  expansion 

L 

denotes 

generalized  Laguerre  expansion 

RC 

denotes 

recursively  via  cumulants 

DM 

denotes 

directly  via  moments 

RM 

denotes 

recursively  via  moments 

Table  F-l.  Shorthand  Notation 

Then,  for  example,  the  combination  PHRC  means  that  this  program  yields  the 
cumulative  or  exceedance  distribution  function  in  terms  of  a  Hermite 
expansion,  by  means  of  expansion  coefficients  determined  recursively  via 
cumulants.  The  eight  programs  listed  here  are,  in  order, 
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PHRC 

Figures  7  and  8 

pHRC 

Figures  7,  9,  and  10 

PHDMandRM 

Figures  7  (and  8) 

pHDMandRM 

Figures  7  (and  9,  10) 

PLRC 

Figures  19  and  20 

pLRC 

Figures  19,  21,  and  22 

PLDMandRM 

Figures  19  (and  20) 

pLDMandRM 

Figures  19  (and  21,  22) 

Table  F-2.  Program  Abbreviations 

The  combination  DMandRM  means  that  this  program  gives  the  expansion 
coefficients  directly  via  moments  as  well  as  recursively  via  moments;  the  user 
must  select  the  procedure  of  interest. 

The  only  input  statistics  we  have  given  a  listing  for  here  is  the  shot 
noise  process  used  in  examples  D  and  I;  in  particular,  the  cumulant  and  moment 
routines  are  listed  at  the  very  end  of  PHRC  and  PHDMandRM,  respectively.  The 
figure  references  given  in  table  F-2  indicate  where  each  particular  program 
was  used  in  this  report;  the  parenthetical  references  are  alternative  ways  of 
generating  those  figures.  The  remaining  figures  in  this  report  require  that 
the  cumulant  and  moment  subroutines  be  replaced  by  the  appropriate  statistics 
of  interest. 
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To  save  space,  no  subroutines  are  listed  more  than  once;  instead, 
comments  are  made  indicating  where  the  needed  routines  are  located,  according 
to  the  coding  in  table  F-2.  For  example,  in  program  PHDMandRM,  function 
subprogram  FNPhi,  line  570,  the  comment  is  made  that  this  routine  has  already 
been  listed  in  PHRC. 

We  now  explain  some  of  the  details  of  the  PHRC  program,  as  an  example,  so 
that  a  user  can  apply  these  techniques  and  routines  to  his  particular 
problem.  The  user  must  specify  M  in  line  30,  which  is  the  maximum  order  of 
approximation  desired,  or  the  number  of  cumulants  or  moments  that  can  be 
calculated.  The  notation  DOUBLE  in  line  40  denotes  INTEGER  variables.  The 
user  must  select  a  and  e  in  lines  130,140;  if  they  are  chosen  equal  to 
ciq,Sq  which  have  been  computed  in  lines  110,120,  then  expansion 
coefficients  a^  =  ^  or  equivalently  b^  =  \>2  ~  0*  However,  this 
choice  is  recommended  only  as  a  starter  on  the  search  in  a,  b  space. 

The  CALL  in  line  150  is  to  the  subroutine  which  calculates  the  expansion 
coefficients  for  a  Hermite  series,  recursively  via  cumulants,  as  can  be 
deciphered  from  the  abbreviated  subroutine  title.  The  expansion  coefficients 
{bn*\are  calculated  and  the  running  sum  of  b^  is  calculated,  both  of 
which  are  printed  on  the  CRT  vs  n.  Also,  a  plot  of  the  expansion  coefficients 
jbn\is  made  in  this  subroutine,  from  which  the  user  must  decide  on  the  order, 

N,  to  employ  in  the  approximate  cumulative  and  exceedance  distribution 
function;  alternatively,  he  can  reject  the  sequence  of[bn}so  obtained,  and 
re-run  the  program  with  different  a,B  values. 
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When  a  satisfactory  o,B  pair  is  obtained,  the  limits  on  the 

range  of  arguments  of  the  distribution  must  also  be  specified;  this  selection 
is  aided  by  the  print-out  of  the  center  and  rms  width  of  the  density  under 
investigation.  A  plot  of  100  values  of  the  cumulative  and  exceedance 
distribution  functions  is  then  made  on  a  logarithmic  ordinate.  The  various 
subroutines  are  self-explanatory  and  are  keyed  to  the  equation  numbers  in  this 
report. 
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PROGRAM  PHRC 


10 
20 
3G 
4  0 
50 
60 
7G 
80 
90 
1  0  0 
110 
12Q 
130 
140 
150 
160 

1  70 
180 
190 

2  0  0 
210 
220 
230 
240 
250 
260 
270 
280 
290 
3  0  0 
3 1  0 
320 
330 
340 
350 
3  6  0 
370 
330 


'6 TER  PL U S  C O H T  I  N U 0 U S  PART  0 F  S H 0 T  N O  I S E  CDF,  P c  *%  u  )  T R  7 3 77,  FIG U R E  8 

COEFFICIENTS  OF  HERMITE  EXPANSION  FOUND  RECURSIVELY  VIA  CUMULANTS 

M=80  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER  OF  CUMULANTS  REQUIRED 

DOUBLE  M , I , N , K  !  INTEGERS  <  2A31  =  2,147,483,648 

RED  I  M  Cum  <  0  :  M  >  ,  A  <  0  :  M  >  ,  He  <  0  :  M  } 

REAL  C  u  fi'i  C  0  :  1  0  0  >  ,  A  <  0 :  1  0  8  >  ,  H  e  <  8  :  1  0  0  >  ,  p  <  0  :  1  0  0  ) 

P0  IS  STEP  AT  ORIGIN 


CALL  C  u  m u  1  an  t  s  <  M ,  P  0 ,  C  u  rn  <  *  )  > 

C  e  n  t-  e  r  =  C  u  m  <  1  > 

R  2  =  C  u  m  ( 2  > 

R  rri  s  =  S  Q  R  (.  R  2  ) 

A  1 p h a0  =  C e n t er 
B  e  t  -=l0  =  R  m  s 

H  1  p  h  a = C  e  n  i  e  r 
B  e  t  a=  R  rn  s  *  1 . 5 

CALL  C  o  e  f  f  h  r  i  a_c  u  m  (  fl ,  A  1  p  h  a ,  B  e  t  a ,  C  u  m  (  #  >  ,  A  <  *  }  ) 
P  R  I  N  T  "  C  e  n  t  e  r  =  M  ;  C  e  n  t  e  r 
PRINT  "  Pros  =  “  ;  R  rn  $ 

F 1  =  1 . /SQR <  2 • *P I > 

INPUT  "ORDER  AND  LIMITS: " , N, U1 , U2 
PRINT  "ORDER  AND  LIMITS: " , N; U1 ; U2 
D u  =  < U  2 - U 1 >/100. 

PLOTTER  IS  "GRAPHICS" 

GRAPHICS  ON 
WINDOW  U1 , U2, -10. , O. 


CENTER  OF  PDF  pcCu) 

MEAN  SQUARE  SPREAD  OF  pc  <u) 
R  M  S  S  P  R  E  A  D  0  F  p  c  <  u  ) 

THE  C  H  0  I C E  S  A  1 p h a= A  1 p h aO 

B  e  t  a=  B  e  t  aO  W  0  U  L  D  M  A  K  E  A  <  1 


RC 


AND 

=  A  <  2  >  =  0 


FOR  1=0  TO  100 

U=U 1 +Du* I 

T  =  ( U  -  A 1  p  h  a  >  /  B  e  t  a 

CALL  Her mi teCN, T, He<*> > 

S  u  rn  =  0 . 

FOR  K= 1  TO  N 
S  u  m  =  S  u  rn  +  A  <.  K  >  *  H  e  <  K  -  1  > 

NEXT  K 

P  =  A  < 0 >  *FNPh i ( T  > -F 1 *EXP < -  *  5*T*T  > *Sum 
IF  U >  =  0 ■  THEN  P  =  P  +  P0 
P< I ) =  p 

IF  P  >  0  a  THEN  400 
PEN  UP 


PROBABILITY  THAT  RV  <  U 
ADDITION  OF  STEP  AT  ORIGIN 


3  9  0  G 0 T  0  4  1 0 

400  PLOT  U , L G T  < P ) 

410  NEXT  I 

420  PENUP 

4  3 0  F  0  R  1=0  T  0  1 0 0 

440  U=Ul+Du*I 

450  Pl  =  l.-Pc:i> 

460  IF  P 1 > 0 .  THEN  490 

470  PENUP 

4  8  0  G  0  T  0  5  0  0 

490  PLOT  U , L G T < P 1 > 

500  NEXT  i' 

510  PENUP 

520  GOTO  190 

530  END 

540  ! 
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558 

DEF  FNPhiC 

568 

Y  =  H  B  S  <  S  Q  R  ( 

578 

SELECT  Y 

538 

C  H  S  E  <  8 . 

5  9  8 

P= 1 63 1 . 760 

1 . 86435 

39749895425 

6  0  0 

P=3723. 587 

132. 267 

0 1 0 8 3  004974 

618 

>3  =  7542.  479 

13, 0777 

1  0  7  5  0  3  6  2  2 1  6 

6  2  0 

Q  =  372  3 . 507 

1349.3465612844574 

630 

Phi =.5*EXP 

648 

L-  ft  y  E  \  2  6 . 6 

6  5  0 

P=2 . 973865 

0 1984 9 

72678426746 

6  6  0 

Q=3 . 369875 

: . 8  4395 

19273551290 

678 

Phi =. 5+ EXP 

630 

CASE  ELSE 

6  9  0 

P  h  i  =  0 . 

7  0  0 

END  SELECT 

710 

IF  X>0.  THI 

728 

RETURN  Phi 

7  3  0 

FNEND 

748 

! 

758 

S  U  E  H  e  r  rn  i  t  > 

760 

DOUBLE  K 

770 

H  e  <  0  )  =  1  . 

738 

He-;:  i  >=:«: 

798 

FOR  K=2  TO 

8  0  0 

He  <  K  >  =X*He 

3  1  8 

NEXT  K 

3  2  0 

SUBEND 

330 

! 

840 

S  U  B  M  o  rn  n  t  1 

358 

DOUBLE  K , J 

8  6  0 

REAL  Mom0 

370 

Morn  <  8  >=  Mo  mi 

8  3  0 

FOR  K = 1  TO 

890 

T=  1  . 

900 

S=Curn  (  K 

9  1  0 

FOR  J=i  TO 

9  2  G 

T  =  T  *  <  K  -  J  >  / , 

930 

S=S  +  T  *  C  u  rn  <  1 

9  4  0 

NEXT  J 

950 

Morn  <  K  )  =S 

960 

NEXT  K 

970 

SUE END 

938 

! 

PROGRAM  PHRC  (cont'd) 

!  HfiRT ,  page  146,  #5783  &  #572= 


eq.  41 


•263 1 +  Y*< 86 . 8327622 1 1 943595 1 +  Y* 
14937+ Y *<6753. 2169641 1848589+Y* 
•  0  y  7  2  +  Y  *  Or!  1  7 , 622386304544077  +  Y  *£ 
■4G55+Y*  < 1 5802 . 5359994028425  +  Y* 


3 9 6 03481623505415+ Y *  <  2 . 26852852876732697  +  Y >  >  > ) > 


H  e  n  (  X  ) 


eq. 


eq.  ft - 6 
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PROGRAM  PHRC  (cont'd) 


9  9  0  SUE  Cum  ri  t  _o  i  a_m  o  m  n  t  <  D  0  U  B  L  E 

I860  DOUBLE  K , J 

1010  RERL  Mom0 

1020  ■Mom0  =  f1om<0>  • 

1 0  3  0  C  u  m  <  0  >  =  L  0  G  (  M  o  m  @  > 

1040  FOR  K= 1  TO  M 

1850  T=  1  . 

1 0  €  0  3 = M  o  m  £  K  ) 

1070  FOR  J=1  TO  K - 1 

1080  T=T*<K-J)/J 

1090  S  =  S  -  T  *  M  o  m  ( J  >  *  C  u  rn  <  K  -  J  > 

1100  NEXT  J 

1110  C  u  m  < K  >  =  S  M  o  m  O 

1120  NEXT  K 

1130  SUBEND 

1140  ! 

1150  s U B  Coe f f h r m i a c u m ( D 0 U B L E 

1160  ALLOCATE  BC0:M> 

1170  DOUBLE  K , J , Mx 

1130  F  =  E  e  t  a*  Bet.  a 

119  @  C  u  m  <  1 )  = C  u  m  <  1 }  -  fl  1  p  h  a )  ✓  B  e  t  a 

1 2  0  0  C  u  m  <  2  )  =  C  u  rn  <  2  )  /  F  - 1 . 

1210  FOR  K  =  3  TO  M 

1220  F=F*Beta*<K-l> 

123  0  C  i  j  rn  ( K )  =  C  u  rn  (  K  >  /  F 

1248  NEXT  K 

12  50  fl  (  0  )  =  B  (  0  )  =  E  X  P  (  C  u  rn  (.  0  )  > 

1260  F= 1 . 

1270  FOR  K= 1  TO  N 

1280  S = 0 . 

1290  FOR  J=1  TO  K 

1 3  0  0  S  —  3  +  C  u  rn  r.  J ^  H  r,  K  —  J  ) 

1310  NEXT  J 

1320  flCK)=S/K 

1330  F=F*K 

1340  B<K)=fl<K)*SQR<F> 

1350  NEXT  K 

1360  Mx=Mx+10 

1370  IF  M x < M  THEN  1360 

1380  Threshold*-?. 

1390  T2=Thresho 1 d*2 . 

1400  V=10.  •"'Threshol  d 

1410  G I N  I T 

1420  PLOTTER  IS  "GRAPHICS" 

1430  GRAPHICS  ON 

1440  WINDOW  0. , FLTCMx> , T2, 0. 

1450  LINE  TYPES 


M, PEAL 


M, REAL 


i 

i 


M o rn <*)  ,  Cum(*  )  ) 


fl 1 pha, Bet  a, Cum <  *  > ,  H ( 


MODIFIED  NORMALIZED 
CUMULANTS  FOR  K=1 

N 0 R MALI2ED  C U M U L fl N T 


eq.  63 

q .  6  2 
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PROGRAM  PHRC  (cont'd) 


1460 

FOR  J — 0  TO  Mx  STEP  10 

1470 

MOVE  J , T 2 

1480 

DRAW  J  j  0 . 

1490 

NEXT  J* 

1500 

FOR  J  =  T 2  TO  0 

1510 

MOVE  0,,J 

1528 

DRAW  M  x , J 

1530 

NEXT  J 

1540 

PENUP 

1550 

LINE  TYPE  1 

1560 

IMAGE  4D ,  2 < 4X , M .  1 7DE  > 

1570 

PRINT  n  K 

1580 

S  u  rn = 0  . 

1590 

FOR  K  =  O  TO  M 

1  6  G  0 

B  =  BOO 

1610 

S  L4  ru  =  S  u  ni  +  B  B 

1620 

PRINT  USING  1560; K, B, 

1  630 

IF  B < V  THEN  1660 

1  640 

Y  =  L  G  T  (  B  ) 

1650 

GOTO  1700 

166  0 

IF  E  > - V  THEN  1690 

1670 

Y  =  T  2  -  L  G  T  <  -  B  > 

1  680 

GOTO  1700 

1  690 

Y  =  T  h  r  t-  s  h  o  1  d 

1  7  0  0 

PLOT  K , Y 

1710 

NEXT  K 

1720 

PENUP 

1730 

SUBEND 

1740 

j 

1750 

S  U  B  C  u  m  u  1  an  t  s  <  D  0  U  E  L  E 

1760 

Over  1 ap  =  6 , 2 

1770 

S i g m aa= 1 , 

1  780 

P  0  =  E  X  P  (  “  0  'v*  *z  r  1  ap  ) 

1  790 

ALLOCATE  Mom<0: M> 

1  8  0  0 

DOUBLE  K 

1 8 1 0 

S  =  S  i  g  rn  aa*  S  i  g  rn  a  a 

1820 

Cum < 0  > =0 . 

1  8  3  0 

C  u  rn  < 1  >  =  0  e  r  1  ap  *  S  i  g  rn  a  a 

1 8  4  0 

C  u  rn  ( 2 )  =  0  m  •=•  r  1  ap  *  S  *  4 .  x  3 

1850 

FOR  K  =  3  TO  M 

1  8  6  0 

C  u  rn  (  K  >  =  C  u  rn  <  K  -  2  >  *  S  *  K  *  K 

1  870 

NEXT  K 

1  8  8  0 

CALL  M  o  rn  n  t  v  i  a  ■:  *.4  rn  n  t  ( 

1890 

M  o  rn  <  G  >  =  M  o  rn  <  0  >  -  P  0 

1  9  0  0 

CALL  C  u  rn  n  t  u  i  a  m  o  rn  n  t.  < 

1910 

SUB END 

BOO 


Sum 


ii 


S  u  rr» 


M 5  REAL  P@,Cum<*>>  !  SHOT  NOISE  eqs.  147-158 

!  Fi  V  .  H  0  .  P IJ  L  S  E  S  S  EC  *  fl  V  ERflGE  P  U  L  S  E  D  U  RfiTI  0  N 
!  PARAMETER  OF  RAYLEIGH  AMPLITUDE  PDF 
!  P R 0 B  A E  I  L  I  T  Y  OF  Z E R 0  fl M P L I T U D E  0 F  S H 0 T  N 0  I S E 
!‘  ARRAY  FOR  MOMENTS 


*.25*PI*SQR< . 5  *  P I ) 


/  (.  K  +  1  > 


M  ,  C  u  m  (  *  >  j  M  o  rn  <  *  )  ) 

!  M 0 M E N T  C 0 R R E C T  I  0 N  F 0 R  I  M P IJ L S E  AT  0 R  I  G I  N 

M  ,  M  o  m  <  *  >  ,  C  u  m  (  *  )  ) 
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PROGRAM  pHRC 


10  ! 

CONTINUOUS  PART  OF  SHOT  NOISE  PDF, 

p  c  <  u  ) 

28  ! 

COEFFICIENTS  OF  HERMITE  EXPANSION 

FOUND  RECURSIVELY  V 

30 

M  =  8 0  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER 

OF  Cl 

40 

DOUBLE  M  ?  I , N  ?  K  '  ! 

INTEGERS  <  2 A 

31  =  ; 

50 

RED  I M  Cum < 0 : M  > , R < 0 : M > , He < 0 : M > 

60 

REHL  Cum  < 0 :  1 80  > ,fl(8:  1 00 > , He < 0 :  1 80 > 

70 

CALL  C  u  rn  u  1  an  t  s  (.  M ,  P  8 ,  C  u  m  <  *  >  >  ! 

P0  IS  STEP  AT 

0  R I  G 

8  0 

C  e  n  t  e  r  =  C  u  rn  C  i  >  | 

CENTER  OF  PDF 

p  c  <  u 

90 

R  2  =  C  u  rn  <  2  >  ! 

MEAN  SQUARE  S 

PREflD 

100 

Rrns  =  SQRCR2>  ! 

RMS  SPREAD  OF 

p  c  C  u 

110 

R  1  p  h  a  8  =  C  e  n  t  e  r  j 

THE  CHOICES 

A  1  p  h  ■ 

120 

Bet  a0  =  Rms  j 

B  e  t  a=  B  e  t  a  O 

WOULD 

130 

R  1  p  h  a= C  e  n  t  e  r 

140 

B  e  t  a=  R  rn  s  #  1 . 5 

150 

CALL  C  o  e  f  f  h  r_v  i  a_c  u  rn  C  M ,  R 1  p  h  a ,  B  e  t  a , 

C  L-i  rn  )  ,  AC*)  ) 

1 

168 

P R  I M T  11  C e n  t  er  =  " ;  C e n t  e r 

170 

PRINT  "  R  rn  s  =  "  ;Rms' 

1 8  0 

F  1  =  1  .  /  (.  B  e  t  a*  S  Q  R  (  2  •  * P  I  )  ) 

190 

INPUT  "ORDER  AND  LIMITS:  ",  N  ,  1J 1 ,  IJ 

2 

2  0  0 

PR  I  NT  "  ORDER  AND  L  I M  ITS:",  N ;  U  i  ;  IJ2 

210 

D  u  =  <  IJ  2  -  U  1  >  ••••'  1  0  0  . 

220 

PLOTTER  IS  "GRAPHICS" 

230 

GRAPH  I C  S  0 N 

240 

WINDOW  U 1 , U 2 , O .  , ,  15 

250 

G R ID  6. , . 0 3 

2  6  0 

FOR  1=0  TO  100 

270 

1 J  =  IJ  1  +  D  L4  *  I 

280 

T  - IJ  -  A  1  p  h  a  )  •••'■  B  *=■  t  a 

290 

CALL  Her-rn  i  t  e  <  N  ,  T  ,  He <  *  }  > 

3  0  0 

S  u  rn  =  A  (  8  ) 

3 1  0 

FOR  K= 1  TO  N 

3  2  0 

S  u  rn  —  s  u  rn  +  A  (.  K )  H  e  ( K  ) 

330 

NEXT  K 

340 

P  =  F  1  *  E  X  P  (  —  .  5  *  T  *  T  )  *  S  u  rn  | 

PDF  OF  RV  AT  IJ 

350 

PLOT  U , P 

360 

NEXT  I* 

370 

PEN  UP 

380 

GOTO  190 

390 

END 

400 

!  USE  ROUTINES  IN  PHRC 

TR 


FIGURE  9 


AND 


RC 
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PROGRAM  PHDMandRM 

10  !  STEP  PLUS  CONTINUOUS  PORT  OF  SHOT  NOISE  CDF,  P < u >  5  COEFFICIENTS  OF 
28  !  HERMITE  EXPANSION  FOUND  DIRECTLY  VIA  MOMENTS  OR  RECURSIVELY  VIA  MOMENTS 
30  M  =  S 8  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER  OF  MOMENTS  REQUIRED 

40  DOUBLE  M,I,N,K  !  INTEGERS  <  2 -'31  =  2,147,483,643 

50  RED  I M  Mom < O : M  > , A < 0 : M > , He  < 0 : M > 

6  U  REAL  M  o m ( 0 :  1 0 0 ) , A ( 0 :  1 0 0 ) , H  e ( 0 :  1 0 0 > ?  p ( 0 :  1 0 0  ) 

70  CALL  Moment s < M , P0 , Mom <  *  )  )  !  P0  IS  STEP  AT  ORIGIN 

8G  C e n t e r  =  M o m < 1 ) / M 0 m < 0 >  1  CENTER  OF  PDF  pc<u> 

90  R2  =  Mom  (.  2  >  /Mom  <  0  >  -Cent  er*Cent  er  !  MEAN  SQUARE  SPREAD  OF  pc<u> 

100  Rms  =  SQR  < R2  >  !  RMS  SPREAD  OF  p c<u) 

1  1  U  A  1  p h a 0 =  L e n t  er  1  T H E  C H 0 1 C E S  A  1  p h a=  A  1  p h a0  A N D 

120  Bet a0=Rms  !  Beta=Beta8  WOULD  MAKE  A < 1 >  =  A < 2 > =0 

1 3  U  A 1  p  h  a  -  C  e  n  t  e  r 

140  B et a=Rms* 1 . 5 

1  5  0  CALL  C  o  e  f  f  h  d_o  i  a_m  o  m  <  M  ,  A  1  p  h  a ,  B  e  t  a ,  M  o  m  •:!  * > ,  A  <  *  >  >  !  D  M 

1 6  U  !  CALL  C o e f  f  h r u i a_m  o m ( M , A 1 p h a , B e  t a , M o m ( *  >  ,  A ( *  > >  !  R M 

1 7  0  P  PINT  " C e  n  t e r  =  “ ; C  e n t  e r 

180  PRINT  " Rms  = " 5 Rrns 

190  F 1  =  1 . / S Q R  < 2 . *P  I  > 

200  INPUT  "ORDER  AND  L I M I TS : " , N , U 1 , U2 

2 1 0  PR  I  NT  " ORDER  AND  L I M ITS s " , N ; U 1 ; U2 

2  2  0  D  u  =  <  I  J  2  -  U 1  )  •••- 1  0  0  . 

230  PLOTTER  IS  "GRAPHICS" 

2  4  0  G  R  A  P  H I C  S  0  N 

250  WINDOW  U1,U2,-10. ,0. 

2  6  0  G R I D  D u *10. ,  1 . 

270  FOR  1=0  TO  100 

280  U=Ul+Du*I 

2  9  0  T  =  (  U  -  A 1  p  h  a  )  •••''  B  e  t  a 

300  CALL  Her mi te<N, T, He<*> ) 

3 1 0  S  u  m  =  0  . 

320  FOR  K  = 1  TO  N 

3  3  0  S  u  rn  =  S  u  m + A  ( K  >  *  H  e  '■  K  - 1 ) 

340  NEXT  K 

350  P  =  A < 0 > *FNPh i < T) -F 1 *EXP ( - . 5*T *T)*Sum  !  PROBABILITY  THAT  RV  <  U 

360  IF  U>=0.  THEN  F'=P  +  P0  !  ADDITION  OF  STEP  AT  ORIGIN 

370  PCI)=P 

380  IF  P > 0 •  THEN  410 

390  PEN UP 

400  GOTO  420 

410  PLOT  U , LGT ( P > 

420  NEXT  I 

430  PENUP 

440  FOR  1=0  TO  100 

450  U  =  U 1 +  D u *  I 

460  P 1=1. -PCI) 

470  IF  P 1 > 0 .  THEN  500 

480  PENUP 

490  GOTO  510 

50O  PLOT  U, LGT CPI) 

510  NEXT  I 

520  PENUP 

5  3  0  G 1 J  T  0  2  0  0 

540  END 

550  ! 
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560 

DEF  FNPhiQO 

570 

!  LISTED  IN  F'HRC 

740 

FHEHD 

750 

i 

760 

SUE  Hermi te (DOUBLE  N 

770 

!  LISTED  IN  PHRC 

830 

SUBEND 

.  8  4  0 

j 

"  85  9 

SUE  Hermit e  i (DOUBLE 

8  6  0 

DOUBLE  K 

870 

Hi  < 0 >  =  1  . 

8  8  0 

Hi(l >  =  X 

390 

FOR  K  =  2  TO  N 

9  0  0 

Hi  00=X*Hi  <K-1  >  +  (K-l 

910 

NEXT  K 

920 

SUBEND 

9  3  0 

1  1 

9  4  0 

S  U  B  M  o  m  n  t  u  i  a  c  u  m  n  t, 

950 

!  LISTED  IN  PHRC 

1070 

SUBEND 

1  080 

j 

1  090 

S  U  B  C  o  e  f  f  h  d  m  i  a  m  o  m  < 

1  100 

ALLOCATE  He<0: M > , F ( O 

1110 

DOUBLE  K , J ,  f  1  x 

1  120 

CALL  Hermi teC N, -Alp h 

1130 

T  =  F  (  0  >  =  1  , 

1  140 

F  0  R  K  =  1  T  0  M 

1  150 

F=F(K)=F(K-1 )  *K 

1160 

T  —  T  *  B  e  t  a 

1  170 

H  e  <  K  >  =  H  e  C  K  >  F 

1180 

M  o  m  <  K  >  =  M  o  rri  (  K )  /  <  F  *  T ) 

1190 

NEXT  K 

1  2  0  0 

FOR  K = 0  TO  M 

1210 

3  =  8 . 

1220 

FOR  J=0  TO  K 

1230 

y  =  S  +  H  e  <.  J  ?  $  M  o  m  (  K  —  J  ) 

1240 

NEXT  J 

1250 

A  <  K  >  =  S 

1260 

NEXT  K 

1270 

MAT  F=SQR< F ) 

1280 

MAT  E  =  ft . F 

1  290 

M  x  =  M  x  +  1  8 

1  3  0  0 

IF  Mx< M  THEN  1290 

1  3  1  0 

T h r  e •=■  ho  1  d  =  - 7  . 

1320 

T  2  =  T  h  r  e  s  h  o  1  d  *  2  . 

1330 

V=10. AThresho 1 d 

1340 

G I  N I  T 

1350 

PLOTTER  IS  "GRAPHICS 

1 3  6  0 

GRAPHICS  ON 

1370 

W I NDOW  0 .  , FLT ( Mx > , T2 , 

1  380 

LINE  TVPE3 

PROGRAM  PHDMandRM  (cont'd) 

!  HART,  page  140,  #5788  &  #5725 

,  R E A L  X,  )  !  H e  f  n ^ X 


!  Hi  >'n<X)  =  <-i  >An  He--n <  i  >0  eq.74- 
!  MODIFIED  HERMITE  POLYNOMIAL 


NORMALIZED  HERMITE  POLYNOMIALS;  eq.  68 
NORMALIZED  MOMENTS  re  Beta;  eq.  69 
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PROGRAM  PHDMandRM  (cont'd) 


1390  FOR  J=0  TO  Mx  STEP  10 

1400  MOVE  J,T2 

1410  DRAW  J,0. 

1420  NEXT  j’ 

1430  FOR  J=T2  TO  0 

1440  MOVE  0.,J 

1450  DRAW  Mx , J 

1460  NEXT  J 

1470  PEHUP 

1480  LINE  TYPE  1 

1490  IMAGE  4D, 2<4X, M. 17DE) 

ISO©  PRINT  "  K  '  BOO  Sum 

1510  Sum=0 . 

1520  FOR  K=0  TO  M 

1530  B  =  BOO 

1540  s  u m = S u m  +  B  * B 

1550  PRINT  USING  1490; K,B, Sum 

156©  IF  B < V  THEN  1590 

1570  Y  =  LGT ( B  > 

1580  GOTO  1630 

1590  IF  B > - V  THEN  1620 

1680  Y=T2-LGT  (. -B  ) 

1610  GOTO  1630 

1620  Y  —  T  h  r  t  '=■  h  o  1  d 

1630  PLOT  K , Y 

1640  NEXT  K 

1650  PENUP 

1660  SUBEND 

1670  ! 

1680  S  U  B  C  o  e  f  t"  h  r_y  i  a_m  o  m  <  D  0  U  B  L  E  M ,  REAL  A 1  p  h  a ,  E  e  t  a ,  M  o  m  <  *  > ,  A  <  *  > 

1 690  ALLOCATE  Hi < 0: M ) , F  C 0 : M  > , B  < 0 : M ) 

1700  DOUBLE  K , J , M x 

171©  CALL  Herm i t e_i <M , A  1 ph a/ Beta, Hi ( * ) > 

1720  T  =  FOO>  =  1. 

1730  FOR  K = 1  TO  M 

1740  F=F<K>»F<K-1)*K 

1750  T=T*Beta 

1760  Hi < K >  =  H i < K > / F  !  NORMALIZED  MODIFIED  HERMITE  POLYNOMIALS 
1770  Mom<K)=MomOO/<F*T>  !  NORMALIZED  MOMENTS  re  Bet 

1780  NEXT  K 

1790  FOR  K = 0  TO  M 

1 8  0  0  S = M  o  m  (  K  ) 

1810  FOR  J=1  TO  K 

1820  S=S-Hi < 

1830  NEXT  J 

1840  AOO=S 

1850  NEXT  K 

i860  MAT  F  =  S  Q  R ( F  ) 

1 8  r  0  MAT  B  =  ft  *  F 
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PROGRAM  PHDMandRM  (cont'd) 


1880 

M  x = M  x  +'10 

1890 

IF  Mx< M  THEN  188.0 

1  9  0  0 

Threshol d=-7. 

1910 

T2  =  Thrs-s-ho1  rf*2. 

192  0 

V  =10.  ''Threshold 

1930 

G  I  N  I  T 

1940 

PLOTTER?  IS  "GRAPHICS" 

1950 

GRAPHICS  ON 

196  0 

Window  0. ,  fltcmx)  ,T2, 

0  . 

1970 

LI'NE  TYPE  3 

1980 

FOR  J=0  TO  Mx  STEP  10 

199  0 

MOVE  J ,  T  2 

2  0  0  0 

DRAW  J,0. 

2010 

NEXT  J 

2020 

F  O  R  J  =  T  2  T  O  0 

2  0  3  0 

MOVE  0.,J 

2  0  4  0 

BRAN  Mx, J 

f 

2050 

NEXT  J 

2  0  6  0 

PEN  UP 

2  9  7  0 

LINE  TYPE  1 

2  0  8  0 

IMAGE  4D, 2<4X, M. 17DE) 

2  0  9  0 

PRINT  H  K 

B 

<  K  > 

2  1  0  0 

S  u  rn  =  0  . 

2110 

FOR  K  =  8  TO  M 

2120 

B  =  E  <  K  > 

2130 

Sum=Sum+B*B 

2140 

PRINT  USING  2080; K,Bf 

S  u  rn 

2150 

IF  E < V  THEN  2130 

2160 

Y  =  LGT  < B  > 

2170 

GOTO  2220 

2 1  8  0 

IF  E > - V  THEN  2210 

2190 

Y  =  T  2  -  L  G  T  <  -  B  > 

2200 

GOTO  2220 

2210 

Y  —  T  h  r  e  s  h  o  1  d 

222  0 

PLOT  K,Y 

2230 

NEXT  K 

2240 

PEN  UP 

2250 

SUB END 

2268 

! 

2  2  7  0 

SUB  Moment  s  < DOUBLE  M, 

REAL 

P  0  ,  C  u  rn  <  *  >  ) 

2  2  8  0 

Oyer  1 ap  =  6 . 2 

i 

AV.  NO.  PULS 

229  8 

S i g m aa= 1 , 

I 

PARAMETER  OF 

2  3  0  0 

P  0  —  E  X  P  '•!  -  0  m  e  r  1  ap  ) 

1 

PROBABILITY 

2310 

ALLOCATE  CumOIi:  M> 

i 

A  R  R  A  Y  F  0  R  C  U 

2320 

DOUBLE  K 

2330 

o  =  o  i  q  fn  a.  =L*  S  i  g  m  =l  =l 

2340 

C  u  m  <  0  >  =  0  . 

2350 

C  u  m  <  1  )  =  0 e  r  1  ap  *  S  i  g  m  a  a 

* .  25* 

P I  *SQR  (. .  5 * P I 

2360 

C  u  rn  <  2  >  =  0  m  e  r  1  a  p  *  S  *  4 .  /  3 

« 

2370 

FOR  K  =  3  TO  M 

2380 

C  u  m  <  K  )  =  C  u  rn  <  K  -  2  >  *  S  *  K  *  K 

x  <  K  + 1 

it 

2390 

NEXT  K 

2  4  0  0 

CALL  M  o  rn  n  t  v  i  a  c  u  rn  n  t  < 

M  ?  C  u  rn 

<  *  >  ,  M  o  rn  <  *  >  > 

24  10 

M  o  rn  i  0 )  =  M  o  rn  C  0  )  -  P  0 

i 

MOMENT  CORRE 

2426 

SUE END 

:•  u  m 


!  SHOT  NOISE 


147-150 
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PROGRAM  pHDMandRM 

10  !  CONTINUOUS  PART  OF  SHOT  NOISE  PDF,  pc<u>;  COEFFICIENTS  OF  HERMITE 


20  ! 

EXPANSION  FOUND  DIRECTLY  VIA  MOMENT 

S  OR  RECURS 

IVELY  VIA  MOMENTS 

3G 

M=80  !  MAXIMUM  ORDER  OF  APPRO 

X 

I  MAT  I  ON;  NUMBER  OF  MOMENTS  REQUIRED 

4  0 

DOUBLE  M,I,M,K  * 

INTEGERS  < 

^ 3 1  —  ci,  147,483, 648 

50 

RED  I M  Mom ( 0 : M > ,  A  <  0 : M >  ,  He < 0 : M > 

60 

REHL  Mom <0: 1 @8 > , R < 0 :  1 08 ),He ( 0 : 1 88 > 

70 

U  R  L  L  M  o  rn  e  n  ts(  M  ,  P  0  ,  M  o  rn  (  *  )  )  ! 

P 0  IS  STEP 

AT  ORIGIN 

30 

C  e  n  t  e  r  =  M  o  rn  (  1  )  s  M  o  m  ( 8  )  ! 

CENTER  OF  PDF  pcCu) 

90 

R  2  =  M  o  rn  <  2  >  •••■'  M  o  rn  <  8  >  -  C  e  n  t  e  r  *  C  e  n  t  e  r  ! 

MEAN  SQUARE 

S  P  R  E  A  D  0  F  p  c  (  u  ) 

1  0  0 

R  rn  s  —  S  Q  R R  2 )  ! 

RMS  SPREAD 

0  F  p «:  <  i.4  > 

1  10 

A  1  p  h  a  8  =  C  e  n  t  e  r  ! 

THE  CHOICES 

A  1  p  h  a=  A  1  p  h  a0  A  N  D 

120 

E  ©  t  a8  =  R  rn  s  ! 

B  e  t  a=  E  e  t  au 

WOULD  MAKE  A<1>=A<2)=8 

130 

A  1  p  h  a  =  C  e  n  t  e-  r 

140 

Bet  a=  R rn s *  1 . 5 

150 

CALL  C  o  €•  +’’  f  h  d  y  i  a  rn  o  rn  ( M ,  A  1  h  a ,  E  e  t  a , 

M 

o rn < * )  ,A(*)) 

!  DM 

1  60 

!  U  A  L  L  C  o  ©  +  ’*  f  h  r _ y  i  a  rn  o  rn  {  M  ,  A  1  fi*  h  a ,  B  e  t  a  n 

M 

o  rn  <  *  )  ,  A  (  *  >  > 

!  RM 

170 

P  R  I  N I  "  C  e  n  t  e  r  =  "  ;  C  e  n  t  e  r 

130 

PRINT  "Rms  =u  ;  Rrns 

190 

F  1  =  1  .  /  <  E  e  t  a*  S Q R  <  2  .  *  P I  )  > 

200 

I MPUT  "  ORDER  AND  L  I  M  ITS:",  N  ,  IJ 1  ,  U 

V 

210 

PRINT  "ORDER  AND  LIMITS: " , N; U1 ; U2 

2  20 

D  u  =  <  1 J  2 - U 1 >/100. 

230 

PLOTTER  IS  "GRAPHICS" 

240 

GRAPH  I C S  0 N 

250 

WINDOW  U 1 , U 2 , 8 . ,  . 15 

260 

GRID  6 ■ f  • 0 3 

27G 

F  0  R  1=0  T  0  1 0 8 

230 

U=Ul+Du*I 

290 

T  =  IJ  ~  A  1  p  h  a )  /  B  e  t  a 

3  0  0 

CALL  He r  rn  i  t  e  <  N ,  T  ,  H e  <  *  >  > 

3  1  0 

S  u  rn  =  A  <  0  > 

3  2  0 

FOR  K= 1  TO  N 

330 

S  u  rn = S  u  rn + A  <  K  >  *  H  e  <  K  > 

340 

NEXT  K 

350 

P  =  F  1  *EXP  <  - .  5 * T * T  >  *Surn 

!  PDF 

OF  RV  AT  U 

3  6  0 

PLOT  U,P 

370 

NEXT  I 

330 

PEN  UP 

390 

G  0  T  0  2  8  0 

4  0  0 

END 

410 

!  USE  ROUTINES  IN  PHDM&RM 
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i 

*  , 


C?  * 


PRO-AM  PLf^C 


10  ! 
20  ! 

3  0 

4  0 
50 
60 
7  0 
3  0 
90 

1  0  0 
110 
120 
130 
140 
150 
160 
170 
180 

1  9  0 

2  0  0 
210 
22© 
230 
240 
25© 

2  6  0 
270 
280 

2  9  O 

3  0  0 
3 1  0 
320 
3  3  0 
340 
350 
3  6  0 
370 
3  8  0 

3  9  0 
400 

4  1  0 
420 
430 
440 
450 
460 
470 
4  8  0 
490 
500 
510 
520 
530 
540 
550 
560 


STEP  PLUS  CONTINUOUS  PORT  OF  SHOT  NOISE  CDF,  Pc<u>;  TR  7377,  FIGURE  20 
COEFFICIENTS  OF  GEN.  LfiGUERRE  EXPANSION  FOUftD  RECURSIVELY  VIA  CUMULRNTS 
1*1=70  !  MAXIMUM  ORDER  OF  APPROXIMATION;  NUMBER  OF  CUMULRNTS  REQUIRED 

DOUBLE  M,I,N,K  !  INTEGERS  <  2A31  =  2,147,483,643 

RED  I M  Cum < Q : M ) , A  C  @ : M > , L  < 0 : M  > 


*  P  <  0 :  1  0  0  > 

P 0  I S*  STEP  RT  ORIGIN 
CENTER  OF. PDF  pc  <u) 

MEAN  SQUARE  SPREAD  OF  pcCu) 

R  M  S  S  P  R  E  A  D  0  F  p  c  (  u  ) 

THE  C  H 0  I C  E  S  Hi  p  h a= A  1 p h a 8  A H  D 

Bet a=  B e t a0  WOULD  MAKE  A < 1 >  = A ( 2 > =0 


RC 


REAL  Cum < 0 :  1 00 > , A ( 0 :  1 80 > , L < 0 :  1 00 
CALL  C  u  m  u  1  an  t  s  ( fl ,  P  8 ,  C  u  rn  <  *  >  ) 

C  e  n  t  e  r = C  u  rn  <  1  >  - 

R2  =  Curn<2> 

R  rn  s  =  S  Q  R  <  R  2  > 

A 1  p  h  a0  =  C  e  n  t  e  r  •£  C  e  n  t  e  n  R  2  —  1  ■ 

E  e  t  a  0 = R  2  C  e  n  t  e  r 
A  1  p  h  a=  .  7  4 
Bet a=2 . 1 

CALL  C o e f  f  1  r_M i  a__e u m < M ,  A 1  pi h a ,  B e t  a «,  C u m ( #)  ,  A ( * > 

P  R I N  T  11 C  e  r*i  t  e  r  =  "  ;  C  e  n  t  e  r 
PRINT  "Rffts  =";Rms 
A 1  =  A 1  p  h  a  + 1 . 

01  =  1 . / A  1 

F  1  =  1 .  /  F  N  ij  am  rn  a  v  A  1  ) 

INPUT  " ORDER  AND  LIMITS: " , N, U1 , U2 
PRINT  "ORDER  AND  LIMITS:  "  ,  N  -  U1 ;  IJ2 
Du=<U2-Ul  >  1 0 0 . 

PLOTTER  IS  "GRAPHICS" 

GRAPHICS  ON 
WINDOW  U1 , U2„ -1 1 . , 8. 

GRID  4.  ,  1  ’ 

P (0)= P  0 

PLOT  0.  , L G T < P 0 ) 

FOR  1=1  TO  100 
IJ  =  U1+Di.4*I 
T = IJ Beta 

CALL  L  ag u  e r  r  e <  N - 1 , A  1 , T , L ( *  >  > 

Sum  =  A  < 0  > *FNF 1 < A  1 , T >*6 1 
FOR  K= 1  TO  N 
S  i.4  rn  =  S  u  rn  +  A  <  K  >  *  L  <  K  -  1  >  /  K 
NEXT  K 

P(  I  )  =P  =  P0  +  F  1  * E X P  - T  +  A 1  * L 0 G  <  T  >  >*Sum  !  PROBAB  I  L  I  TV  THAT  RV 

IF  P  >  0 ■  THEN  420 

PEN  UP 

GOTO  430 

PLOT  U , L G T < P ) 

NEXT  I 
PENUP 

F  0  R  I  =  0  T  0  1  8  0 
U=U 1 +Du* I 
P1=1.-P(I) 

IF  P 1 >  0  «  THEN  510 

PENUP 

GOTO  520 

PLOT  U,LGT<P1) 

NEXT  I 
PENUP 
G 0 TO  21 0 
END 

i 
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PROGRAM  PLRC  (cont'd) 


57  U  D  E  F  F N L -=±rii m a i  X  )  !  G arn rn a (  X  )  v  i  a  H fl R T  ,  p aq e  282,  #524 3 

580  DOUBLE  N , K 

590  H  =  I  N T  <!  X  > 

600  R  =  X - N 

610  IF  N > 0  OR  RO0.  THEN  640 

628  PRINT  "FNGammaCX)  IS  HOT  DEFINED  FOR  X  =  " ; X 
630  STOP 

640  IF  R>0.  THEN  670 

6  5  U  G  am  m  a  2  =  1  . 

660  GOTO  740 

678  P  =  4 3 9 . 3 3 0444  8  6  8  0 2  5  6  7  6  +  R 

6  8 0  P  =  8  7  6  2 . 71 U 297852 1 4  8  9  6  +  R 

690  P  =  4  2  3  5  3 . 6895097440896+R 

760  q  =  499 ..  02852662  1 439048-R 

7 1 0  Q =994 0 . 3074 1 508277090-R 

720  0  =  42353. 6  8  9  5  0  9  7  4  4  0  9  0  0  +  R 

7  3  0  G  arn  m  a 2  =  P  /  Q  !  G  am  m  a  <  2  +  R  }  f  o  r  8  <  R  <  1 

748  IF  N >2  THEN  780 

750  IF  N < 2  THEN  830 

7  6  0  G  a  rn  m  a = G  a  m  rn  a  2 

770  RET  U  R  N  G  am  rn  a 

7  8  y  i  j  a  rn  rn  a = G  a  rn  rn  a  2 

790  FOR  K  = 1  TO  N-2 

8  y  y  ij  am  m  a= G  am  rn  a*  ( X  -  K ) 

810  NEXT  K 

820  RETURN  Gamma 

830  R  =  1  . 

840  FOR  K = 0  TO  1-N 

8  5  0  R  =  R*  (  X  +  K  *•' 

868  NEXT  K 

8  7  y  U  a  rn  m  a  =  b  a  rn  rn  a  2  /  R 

8  y  0  RET  U  R  N  G  a  rn  rn  a 

898  FNEND 

9  0  0  ! 

910  DEF  FNF 1 <  R 1 , X  >  !  1 F  1  < 1 ; R 1  +  1 • X ) 

928  DOUBLE  K 

938  T  =  8= 1 . 

940  F 0 R  K =1  TO  2 8 8 

950  T  =  T  *  X  s ( fl 1 + K > 

960  S=S+T 

978  IF  T <  =  1 . £-17*8  THEN  RETURN  S 

988  NEXT  K 

990  PRINT  “200  TERMS  IN  FNF1  fl T " ; 8 1 ; X 

1008  RETURN  S 

1010  FNEND 

1  0  2  0  ! 


. y  .  1 086937529799530  +  8*6 .  74495072459252899 
2  y  0  8  ■  52  {  4  0  1  3  U  7  2  7  9  1  2  +  R  +  F' 

2 y  3 3 6 . 86178926 9 8 8 7 4  +  R P  ) 

189. 498234 1 578280 iS-R*  <  23 . 881551524588125 
1528. 60727377952202+R+Q > 

9 fl  _  Fi fl  q q  —  p  n  s 
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PROGRAM  PLRC  (cont'd) 


1038 

1048 
1  050 
1060 
1  070 
1  030 
1 0  3  0 
1100 
1110 
1128 
1130 
1148 
126  8 
1278 
1  2  3  8 
1298 
1428 
1430 
1448 
1450 
1468 
1470 
1480 
1490 
1500 
1510 
1520 
1538 
1540 
1550 
1560 
1578 
1530 
1  590 
1  6  0  0 
1610 
1  620 
1  6  3  0 
1648 
1658 
1660 
1670 
1  630 
1  6  9  O 
1  700 
1710 
1728 
1730 
1  740 
1750 


<  ‘ 
e  , 

*  i 

TR*  7377  • 

tl 

#  11  « 


*  •  * 

h  U  E  L  ag  u  e  r  r  £•  D  U  U  E  L  E  N ,  R  E  ft  L  Hi  p  h  a  ,  X  ,  L  \  *  )  )  I  e  q  .  3  6 

DOUBLE "K 
R 1  =  R  1  p  h  a-  1  . 

L  <  0  >  =  1  . 

L  t  1  >  =R  1  p h a+  1  .  “ X 
FOR  K  =  2  TO  N 

L  <  K  =  <  <  K+K+fll-X  >  *L<  K- 1  >  -  <  K+R  1  )  *L  C  K-2  >  > /K 
NEXT  K 
SUE END 
! 

S  U  E  M  o  rn  n  t  _  v  i  a__  c  u  rn  n  t  <  D  i  J  U  E  L  E  M ,  R  E  R  L  C  u  m  <  *  > ,  M  o  m  <  *  >  > 

!  LISTED  IN  PHRC 
SUE END 

S U B  C L4 rn n t _ v  i  a_rn o rn n t  ( D 0 U E L E  M,  REAL  M o m  ( *  ) ,  C u m <  * >  > 

!  LISTED  IN  PHRC 
SUE END 

I 

S  U  E  C  o  <=•  f  f  1  i  a_c  u rn  ( D  0  U  E  L  E  M  ,  R  E  fl  L  R 1  p  h  a ,  B  e  t  a ,  C  u  rn  <  *  > ,  R  <  *  >  > 

RLLOCRTE  B < 0 : M > , C <0 : M > ? D < 1 : M > 

DOUBLE  K, J, J1 , Mx 
T  =  B  e  t  a 

Cum  <  1  >  =Curn  <  1  >  /T 
F 0 R  K  =  2  TO  M 
T  =  T  *  B  e t a*  < K  — 1 > 

C u rn < K )  =  C u m  < K ) / T  !  N 0 RMRLIZED  C U M 1J L R N T S  ;  e q  .  6 2 

NEXT  K 
R 1 =R 1 pha+1 . 

FOR  J=1  TO  M. 

J  1  =  J  +  1 
T=  1  . 

S  =  R  1 

FOR  K = 1  TO  J 
T  =  T  *  (  K  -  J  1  )  s  K 
S  =  S  +  T *Cum < K > 

NEXT  K 
D  < J  > =S 
NEXT  J 

R  (  0  >  =  E  (  0  )  —  C  (.  0  )  =  E  X  P  (  C  u  rn  (  0  )  ) 

Q=  1  . 

FOR  K  = 1  TO  M 
S  =  0 . 

FOR  ■  J  = 1  TO  K 
S  =  S  +  D  k  J  )  *  C  k  K  ~  J  ) 

NEXT  J 

C  (  K )  =  C  =  S  K 

Ui  =  U  *  K  •••■' '%  R  1  p  h  a+  K  ) 

fl<K)=C*Q 

E ( K > =C*SQR  < Q  > 

NEXT  K 
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1760 
1770 
1  7  8  0 
1790 
1800 
1810 
1820 
1830 
1840 
l  y  5  0 
I860 
1870 
1  880 
1  890 
1900 
1910 
1920 
1930 
1940 
1950 
1  9  6  0 
1970 
1  980 

1  990 
2000 
2010 
2020 
2030 
2040 

2  0  5  0 
2060 
2070 
2080 
2090 
2100 
2110 
2120 
2130 
2140 
2150 
2  1  60 
2310 


PROGRAM  PLRC  (cont'd) 


Mx=Mx+ 1 0 

IF  Mx<  M  THEN  1760 
T  h  r  e  ■=■  hoi  d  =  -  7 . 

T2=Thresho 1 d*2 . 

Y  =10.  T  h  r  e  s  h  o  1  d 
G I N I T 

PLOTTER  IS  "GRAPHICS " 

GRAPHICS  ON 

WINDOW  0. , FLT <Mx) , T2, 0. 

LINE  TYPE  3 

FOR  J — 0  TO  Mx  STEP  10 

MOVE  J ,  T  2 

DRAW  J,0. 

NEXT  J 

FOR  J  =  T 2  TO  0 
MOVE  0 ■ f J 
DRAW  M x , J 
NEXT  J 
PENUP 

LINE  TYPE  1 

IMAGE  4  D ,  2  <  4 X , M . 17DE> 

PRINT  "  K  *  BOO  Sum" 

S  u  m  =  0  . 

FOR  K  =  8  TO  M 
B = B  (  K  > 

S  u  m  =  S  1.4  m  +  B  *  E 

PR  I  NT  US  I NG  1 960 ; K , E , Sum 
IF  B < V  THEN  2060 

Y  =  L  G  T  <  B  > 

GOTO  2100 

IF  E  > - V  THEN  2090 

Y  =  T  2  -  L  G  T  <  -  E  > 

G 0 TO  21 0 O 

Y  =  T  h  res h  o 1 d 
PLOT  K , Y 
NEXT  K 
PENUP 

SUB END 
! 

SUE  C L4 m u  1  ant  s  <  DOUBLE  M ,  REAL  PO ,  Cum  <  *  >  >  !  SHOT  NO  1 3E 

!  LISTED  IN  PHRC 
SUBEND 
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PROGRAM  pLRC 


1  0 
20 
30 
4  0 
50 
60 
70 
80 
90 
1  0  0 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 


U O N 1 1 N U O U S  P H R T  0 F  S H 0 T  N 0 1 S E  PDF,  p c  < u >  TR  7377,  FI G U RE  21 
COEFFS.  OF  GENERAL.  LAGUERRE  EXPANSION  FOUND  RECURSIVELY  VIA  CUMULANTS 
M  =  ^ y  M  A X  I  M U M  0 R D E R  OF  A P P R 0 X  I  M  A T  I  0 N ;  N U M B E R  0 F  C U M !J L  A N T S  R E Q I J  I  R  E  D 


DOUBLE  M , I , N , K  !  INTEGERS 

RED  I  M  Cum  <  0  :  M  >  ,  A  •:!  0 :  M  > ,  L  <  8  :  M  > 

REAL  C  u  m  0 :  1  O  0  >  ,  A  <  0  :  1  8  0  )  ,  L  <  0 :  1 0  0  > 

CALL  C  L4  m  u  1  an  t  s  C  M ,  P  0 ,  C  u  m  (  # ) 

C  e  n  t  e  r  =  C  u  m  <  1  ) 

R  2  =  C  u  m  <  2  !:■ 

R  m  s  =  S  Q  R  <  R  2  ) 

A  ]  p  h  a  0  =  C  e  n  t  e  r  *  C  e  n  t  e  r  /  R  2  - 1 , 

E  e  t  a@  =  R  2  ■  C  e  n  t  e  r 
A  1 p ha- . 7  4 
Beta® 2. 1 

CALL  C o e  f' f  1  r u  i  a c  u m  {  M ,  A 1  p h a ,  B •=■  t  a ,  C u m  ( *  )  ,  A  (  * ) ) 

PRINT  "Center  =  " ; Center 
PRINT  " R rn s  =";Rms* 

F  1  =  1  .  f  •*.  B  e  t  a*F  N L  am m  a (  A  1  p h  a+  1  .  )  ) 

INPUT  "ORDER  AND  LIMITS: " ,N,U1,U2 
PR  I  NT  "  ORDER  AND  L  I M  ITS:",  N  ;  U  1  :  IJ2 


2^-31  =  2,147,483,648 


P 0  IS  STEP  AT  ORIGIN 
CENTER  OF  PDF  pc<u) 

MEAN  SQUARE  SPREAD  OF  pc<u) 

R  M  S  S  P  R  E  A  D  0  F  p  c  <  u  ) 

THE  CHO I CES  A  1 pha=A 1 pha0  AND 
Bet  a=  B e t a 8  WOULD  MAKE  A  < 1 >  = A  < 2  >  =8 


RC 


210  D  u  =  <  U  2-1.11  >  /  1  0  0  . 

220  PLOTTER  IS  "GRAPHICS" 

230  GRAPHICS  ON 

248  WINDOW  U 1 , U2 , 0 . , . 1 5 

250  GR ID  6 . , . 83 

260  FOR  1=0  TO  100 

270  U=U 1 +Du* I 

280  IF  U<0.  THEN  488 

290  IF  U>0.  THEN  320 

380  PLOT  0 ■ , 0 ■ 

3  1  8  G  0  T  0  4  0  0 

328  T  =  IJ  Beta 

■-*  U  L  A  L  L  L  ag  u  e  r  r  e  (  N ,  A 1  p  h  a ,  T ,  L .( * ) ) 

340  Sum  =  A  <  0  > 

350  FOR  K= 1  TO  N 

360  Surri=Sum  +  A  <  K )  *L  ( K  > 

378  NEXT  K 

388  P  =  Fl*EXPOT+Al  pha*LOG<T>  )*Sum  !  PDF  OF  RV  AT  II 

398  PLOT  U,P 

480  NEXT  I 

418  PENUP 

4  2  8  G  0  T  0  1  9  0 

430  END 

448  !  USE  ROUTINES  IN  PLRC 
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PROGRAM  PLDMandRM 


1  y  !  S T E P  P L U S  C 0 M T I N U 0 IJ S  PART  0 F  S H 0 T  N 0  I  S E  CDF,  P c  < u >  ;  COEFFICIE N T R  0 F 

20  !  GENERALIZED  LAGUERRE  EXPAN.  FOUND  DIRECTLY  AND  RECURSIVELY  VIA  MOMENTS 

! I MAT I ON;  NUMBER  OF  MOMENTS  REQUIRED 
INTEGERS  <  2-'-31  =  2,147,483,648 

, P<6: 10@) 

P0  IS  STEP  AT  ORIGIN 
CENTER  OF  PDF  pc<u> 

MEAN  SQUARE  SPREAD  OF  pc  <u> 

R  M  S  S  P  R  E  A  D  0  F  p  c  (  i.j  ) 

THE  CH 0 ICES  A  1 p h a= fl 1 p h a@ 


3  0 

M  =  70  !  MAX  I 

40 

DOUBLE  M , I , H , K 

50 

R  E  D  I  M  ti  o  m  C  Q :  M  > 

60 

REAL  MomC0: 100 

70 

CALL  Moments  CM 

8  0 

Cent  er=Mom  C  1  '>  / 

90 

R  2 = M  o  m  C  2  >  /  M  o  m  C 

1  00 

Rms=SQRCR2> 

1  10 

fl  1  p  h  a8  =  C  e  n  t  e  r  * 

120 

Bet  aQ=R2/Cent  e 

130 

fll pha=. 74 

140 

B e t  a= 2 .  1 

1  50 

CALL  Coe  ft"  1  d_u 

1  60 

!  CALL  Co  eft"  ]  r  v 

170 

PRINT  “Center 

1  3  0 

PRINT  "Rms  =  “;l 

190 

fl 1 =  fl 1 p  ha+ 1 . 

200 

01=1 . /fl 1 

210 

F 1  =  1 .  •••"  F  N  G  a  m  m  a  C I 

2  2  0 

INPUT  "ORDER 

230 

PRINT  "ORDER  fl! 

240 

D u = < U  2 - U 1 >/100 

250 

PLOTTER  IS  "GRl 

260 

G  R  fl  P  H  I  C  8  0  N 

270 

WINDOW  1-11,112,- 

2  S  0 

GRID  4. , 1 . 

290 

P  (  0  )  =  P  0 

3  0  0 

PLOT  0. , LGTCP0 

3  1  0 

FOR  1=1  TO  100 

3  2  0 

u=ui+Du*r 

3  3  0 

T =1-1/ Bet  a 

340 

CALL  Laguerred 

350 

S  L4  m  =  fl  '•  0  >  *  F  N  F  1  O 

3  6  0 

FOR  K= 1  TO  N 

370 

S  u  rn  =  S  u  m  +  fl  (  K  )  *  L 

330 

NEXT  K 

390 

PC  I )=P  =  P0  +  F1*E: 

4  0  0 

IF  P  >  0  ■  THEN  4: 

410 

PENIJP 

420 

GOTO  440 

430 

PLOT  U,LGTCP) 

440 

NEXT  I 

450 

PENUP 

4  6  0 

F  0  R  1=0  T  0  1 0  0 

470 

IJ  =  U  1  +  Du*  I 

430 

Pl  =  l . -PC  I > 

4  9  0 

IF  F 1 > 0 .  THEN  ! 

5  0  0 

PENUP 

510 

i  j  0  T  IJ  5  •  J  u 

520 

PLOT  IJ,LGTCPn 

530 

NEXT  I 

540 

PENUP 

550 

GOTO  220 

560 

END 

570 

i 

Bet  a=Bet  aO 


,  A  <  *  >  ) 


AND 

WOULD  MAKE  fl(l) =A < 2 ) =0 


DM 

RM 


; Cent  er 


0 . 


L  <  *  )  ) 


PROBABILITY  THAT  R'. 
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PROGRAM  PLDMandRC  (cont'd) 


580 

DEF  F N La arri rn a < 

X  >  !  G  am  rn  a  <  K 

*  y i  a  H  R  R T ,  p aq e  282,  #524 

590 

!  LISTED  IN  PLRC 

9  O  0 

FNEND 

9  1  0 

j 

920 

DEF  FNFltAl, 

y  si 

!  lFl<ljfll+l;X> 

9  3  0 

!  LISTED  IN  PLRC 

1  0  2  0 

FNEND 

1 030 

i 

1  040 

SUE  La guerre 

<  DOUBLE  N , REAL 

R  1  fj  h  a ,  X  ,  L  *  ..i  !  L  n  \  a  1  p  h 

1050 

!  LISTED  IN  PLRC 

1120 

SUE END 

1  1  30 

j 

1140 

S  U  B  M  o  rn  n  t  v  i 

a  c  u  rn  n  t  <  Ii  0  l_l  B  L  E 

M ,  R  E  H  L  C  u  rn  <  *  >  ,  M  o  rn  (  *  >  > 

1150 

!  LISTED  IN  PHRC 

1270 

SUBEND 

1  2  8  0 

! 

1290 

SUE  Coef'fld 

v  i  a  rn  o  rn  <  D  0  U  E  L  E 

M  ,  R  E  R  L  R  1  p  h  a ,  E  e  t  a ,  M  o  rn  <  *  )  , 

1 300 

ALLOCATE  B<0 

:  M  ) 

1  3  1  0 

DOUBLE  K j  K 1 , 

J  ,  fix 

1320 

T=  1  . 

1330 

FOR  K= 1  TO  M 

1 340 

T  =  T  *  c.  A 1  p  h  a+  K 

)  *  B  e  t  a 

!  NORMALIZED  MOMENT 

1350 

M  o  rn  <  K  >  -  M  o  rn  ■:!  K 

>  /  T 

!  Alpha  and  Beta;  e 

1360 

NEXT  K 

1370 

Q=  1  . 

1  3  8  0 

H  <  0  >  =  B  <  0 )  =  M  o  rn  <  0  > 

1  3  9  0 

FOR  K= 1  TO  M 

1400 

K  1  =  K  +  1 

1410 

T=  1 . 

1420 

S  =  M  o  rn  ( 0 ) 

1430 

FOR  J=1  TO  K 

1440 

T=T*< J-Kl >/J 

1450 

S  =  S  +  T  *  M  o  rn  <  J  ) 

1460 

NEXT  J 

1  470 

Q  =  Q  *  • .  R  1  ph  a+  K 

)/K 

1  480 

fl(K)  =  S 

1490 

B  K —  S  #  S  Q  R £  Q 

> 

1500 

NEXT  K 

1510 

Mx=Mx+  1 0 

1520 

IF  Mx< M  THEN 

1510 

1530 

T  h  r  e  s  h  o  1  d  —  -  7 

m 

1540 

T  2 = T  h  r  e  s  h  o  1  d 

*2. 

1550 

V  -lu.  T  h  r  e  s  h  o  1  d 

1560 

G  INI  T 

1570 

PLOTTER  IS  " 

GRAPHICS'1 

1580 

G  R  fl  P  H  I  C  S  0  N 

1590 

WINDOW  O. , FLT < Mx> , T2, 0. 

1  600 

LINE  TYPE  3 
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PROGRAM  PLDMandRM  (cont'd) 


1610 

FOR  J  =  0  TO  Mx  STEP  10 

1  62U 

MOVE  J ,  T  2 

1  630 

DRAW  J,0. 

1640 

NEXT  J 

1650 

FOR  J  =  T 2  TO  O 

1  660 

MOVE  0 . , J 

1  670 

DRAM  Mx, J 

1  680 

NEXT  J 

1  690 

PEN  UP 

1  7  0  0 

LINE  TYPE  1 

1710 

IMAGE  4D, 2C4X, M. 17DE) 

1720 

PRINT  "  K  '  BOO 

1 7  3  0 

S  u  m  =  0  . 

1740 

FOR  K  =  0  TO  M 

1750 

E  =  EOO 

1760 

S  u  rn  =  S  u  m  +  E  *  B 

1770 

P  R I N  T  U S I N G  171 0 ; K , B , S u  m 

1  780 

IF  E  <  V  THEN  1810 

1790 

Y=LGT(E> 

1 800 

GOTO  1850 

1  8  1  0 

IF  B > - V  THEN  1840 

1  820 

Y  =  T  2  -  L  G  T  <  -  B  > 

1  830 

G  0  T  0  1  8  5  0 

1  840 

Y  =  T  h  r  e  s  h  o  1  d 

1  8  5  0 

PLOT  K,Y 

1  360 

NEXT  K 

1870 

PEN  UP 

1  880 

SUBEND 

1  890 

! 

1  9  0  0 

S U B  C o e f  f  1  r _u  i  a_m o iri  <  D 0 U B L E  M,  REAL 

A  1  p  h  a ,  B  e  t-  a ,  M  o  rn '%  *  >  , 

1910 

ALLOCATE  EOH:M) 

1920 

DOUBLE  K , K 1 , J , Mx 

1  930 

T=  1  . 

1  940 

FOR  K = 1  TO  M 

1  950 

T  =  T  fi  1  p  h  a+k  ★  E  *•  t  a  | 

NORMALIZED  MOMENT 

I960 

M  o  rn  (  K  )  =  M  o  rn  C  K  )  s  T  | 

A 1  p  h  a  an  d  E  *=•  t,  a ;  *=• 

1970 

NEXT  K 

1  9  8  0 

Q=  1  . 

1990 

A  <  O  >  =E  <  O ) =Mom  <  O  > 

2  0  O  0 

FOR  K = 1  TO  M 

2010 

K 1 =K+ 1 

2  0  2  0 

T=  1 . 

2  0  3  0 

S  =  Morn,OO-A<0> 

2040 

FOR  J  = 1  TO  K- 1 

2050 

T=T*(J-K1)/J 

2  0  6  0 

S  =  S  -  T  *  A  '■  J  > 

2070 

NEXT  J 

2080 

IF  K  MOD  2=1  THEN  S=-S 

2  0  9  0 

AOO=S 

2  1  0  0 

Q  =  9.  *  ■:!  A  1  p  h  a+ K  )  f  K 

2110 

E < K >  =  S*SQR i Q  > 

2120 

NEXT  K 
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2136 
2140 
2150 
2160 
2170 
2180 
2190 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2280 
2290 
2300 
2310 
2320 
2330 
2340 
2350 
2360 
2370 
2380 
2390 
2400 
2410 
2420 
2430 
2440 
2450 
2460 
2470 
2  4  8  0 
2490 
2500 
2510 
2520 
2530 
2670 


PROGRAM  PLDMandRM  (cont'd) 


Mx=Mx+ 1 0 

IF  Mx<  M  THEN  2130 
T  h  r  e  s  h  o  1  d  =  -  7  , 

T2=Thresho 1 d*2 . 

Y  =10.  T  h  r  e  s  h  o  1  d 
G I  N I T 

PLOTTER  IS  "GRAPHICS11 
GRAPHICS  ON 

WINDOW  0. , FLTCMx) , T2, 0. 

LINE  TYPE  3 

FOR  J=0  TO  Mx  STEP  10 

MOVE  J,T2 

DRAW  J„0. 

NEXT  j’ 

FOR  J=T2  TO  0 
MOVE  0 ■ j  J 
DRAW  Mx, J 
NEXT  J 
PEN  UP 

LINE  TYPE  1 

IMAGE  4 D , 2(4 X ,  M .  1 7 D E > 

PRINT  "  K  BOO  Sum 

S  u  m  =  0 . 

FOR  K  =  0  TO  M 
B  =  B  <  K  > 

S  u  m  =  S  u  m  +  E  *  E 

PR  I  NT  US  I NG  2330 ; K , B , Sum 
IF  B < V  THEN  2430 

Y  =  LGT  < B  > 

GOTO  2470 

IF  E  > - V  THEN  2460 

Y  =  T2-LGTOB> 

GOTO  2470 
Y=Thresho 1 d 
PLOT  K , Y 
NEXT  k’ 

PEN  UP 
SUE END 
! 

SUE  Moment  s < DOUBLE  M , REAL  PO , Mom  < *  > >  !  SHOT  NO  I SE 

!  LISTED  IN  PHDM&RM 
SUB END 
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PROGRAM  pLDMandRM 


18  ! 
20  ! 
30 
40 
50 
60 
70 
80 
90 
1  0  0 
1  10 
120 

1  30 
140 
150 

160  ! 
170 
180 
190 

2  0  0 
210 
2  2  0 
2  3  0 
240 
250 
260 
270 
2  8  0 

2  9  0 
300 

3  1 0 
320 
330 
340 
350 
3  6  0 
370 
380 

3  9  0 

4  0  0 
4  1  0 


CONTINUOUS  PART  OF  SHOT  NOISE  PDF,  pc  0.0;  COEFFICIENTS  OF  GENERALIZED 
LAGUERRE  EXPANSION  FOUND  DIRECTLY  AND  RECURSIVELY  VIA  MOMENTS 
M - 7  U  !  M A  X  I M U M  0  R  D  E  R  0 F  APR  R  0 X  I M  A  T I  0 N ;  N  U M  E E  R  0 F  M  0 M E  N T  S  R E  Q U I R  E  D 


E  e  t  a=  B  e  t  a0 


DOUBLE  M , I , N ?  K 

RED  I M  Morn  C  0  :  M  >  ,  A  <  8 :  M  )  ,  L  <  0  :  M  > 

REAL  Morn  ■:!  0 :  1 00  )  ,  A  <  0 :  1  00  >  ,  L  (  3  :  1  00 
CALL  M  o  rn  e  n  t  s  <  M  ,  P  0 ,  M  o  rn  C  *  >  >  ! 

C  e  n  t  e  r  =  M  o  rn  <  1  >  /  M  o  rn  <  0  >  | 

R  2  =  M  o  rn  (.  2  )  /  M  o  rn  (  0  >  -  C  e  n  t  e  r  *  C  e  n  t  ■=•  r-  ! 

R  rn  s = S  Q  R  ( R  2  )  \ 

A  1  p  h  aU  ~  u  0  n  t-  0-  f"  &  C  ■=■  n  t  e  r*  R  2  —  1  ■  J 

E  e  t  a  0  =  R  2  C  e  n  t  e  r  i 

A  1 pha= . 7  4 
Bet a=2 . 1 

-ALL  i_.  o  eft  1  d  m  i  a _ rri  o  rn  *•.  M ,  A 1  p  h  a ,  B  t  a ,  M  o  rn  #  . j  ?  A  '■!  #  )  ) 

CALL  Coef  f  1  r_M  i  a__rnorn  <  M  ,  A  1  pha.  Bet  a.  Morn  <  *  }  ,  A  <  *  )  ) 
P  PINT  "  C €  n t e r  ~  " ; C e n t e r 
PRINT  11 R rn s  =“  ;Rrns‘ 

F  1  =  1  .  s  <.  E  e  t  a*  F  N  G  am  rn  a  <  A  1  p  h  a+  1  .  >  > 

INPUT  "ORDER  AND  L I M I TS s " , N , U 1 , U2 
PR  I  NT  "  ORDER  AND  L  I  M  ITS  :  "  ,  N ;  I J 1  ;  IJ2 
D  u  =  <  U  2  —  U  1  >  /  1  0  0  . 

PLOTTER  IS  "GRAPHICS" 

GRAPH  I C S  0 N 
WINDOW  U1 , U2, 0. , . 15 
GR ID  6, , . 03 
PLOT  8 • , 0 . 

FOR  1=1  TO  180 
U=U 1 +Du* I 
T  =  U  s  B  e  t  a 


I  N  T  E  L  E  R  S'  2  31  —  2  j,  1 4  7 , 4  8  3 , 64  8 


PO  IS  STEP  AT  ORIGIN 
CENTER  OF  PDF  pc (u) 

MEAN  SQUARE  SPREAD  OF  pc(u) 

R  M  S  S  P  R  E  A  D  0  F  p  c  (  u  ) 

THE  CHO I CES  A 1 pha=A 1 phaQ  AND 

WOULD  MAKE  A < 1 ) = A < 2 > =0 


•DM 

RM 


CALL  L  a  q  u  e  r  r*  e  *•!  N ,  A  1  fi*  ha,  T  ,  L  (  *  ')  > 

S  u  rn  =  A  (  0  ) 

FOR  K=1  TO  N 
S  u  rn  =  S  u  rn  +  A  <  K  >  *  L  <  K  ) 

NEXT  K 

F  =  F  1  * E X P  ( - T  +  fl  1  p h a* L 0 G  <  T  1  >  * S u m  !  PDF  OF  RV  AT  U 

PLOT  U , P 

NEXT  I 

F'EHIJP 

GOTO  200 


END 

USE  ROUTINES  IN  PLDM&RM 


126 


TR  7377 


REFERENCES 


1. 

J.  I.  Marcum,  "A  Statistical  Theory  of  Target  Detection  by  Pulsed  Radar: 

Mathematical  Appendix,  "Research  Memorandum  RM-753,  RAND  Corp.,  Santa 

Monica,  Calif.,  1  July  1948.  Also  in  IRE  Trans,  on  Information  Theory, 

0 

Vol .  IT-6,  No.  2,  April  1960. 

2. 

H.  Cramer,  Mathematical  Methods  Of  Statistics,  Princeton  University 

Press,  Princeton,  N.  J. ,  1961. 

3. 

M.  G.  Kendall  and  A.  Stuart,  The  Advanced  Theory  of  Statistics;  Vol.  1, 

Distribution  Theory,  Hafner  Publishing  Co.,  N.  Y.,  3rd  edition,  1969. 

4. 

A.  H.  Nuttall,  "Under-Ice  Roughness:  Shot  Noise  Model,"  NUSC  Technical 

Memorandum  841208,  31  December  1984. 

5. 

Handbook  of  Mathematical  Functions,  U.  S.  Department  of  Commerce. 

National  Bureau  of  Standards,  Applied  Mathematics  Series  No.  55,  U.  S. 

Government  Printing  Office,  June  1964. 

6. 

4 

A.  H.  Nuttall,  "Recursive  Inter-Relationships  Between  Moments,  Central 

Moments,  and  Cumulants,"  NUSC  Technical  Memorandum  No.  TC-201-71, 

12  October  1971. 

7. 

G.  Szego,  Orthogonal  Polynomials,  American  Mathematical  Society.  Vol.  23. 

Providence,  R.  I.,  3rd  edition,  1967. 

127 


TR  7377 


8.  I.  S.  Gradshteyn  and  I.  M.  Ryzhik,  Table  of  Integrals,  Series,  and 
Products,  Academic  Press,  Inc.,  New  York,  1980. 

9.  A.  H.  Nuttall  and  B.  Dedreux,  "Exact  Operating  Characteristics  for  Linear 
Sum  Of  Envelopes  of  Narrowband  Gaussian  Process  and  Sinewave,"  NUSC 
Technical  Report  7117,  11  January  1984. 

10.  A.  H.  Nuttall,  "Accurate  Efficient  Evaluation  of  Cumulative  or  Exceedance 
Probability  Distributions  Directly  from  Characteristic  Functions,"  NUSC 
Technical  Report  7023,  1  October  1983. 

11.  A.  H.  Nuttall,  "Exact  Performance  of  General  Second-Order  Processors  for 
Gaussian  Inputs,"  NUSC  Technical  Report  7035,  15  October  1983. 


128 


TD  7377 


INITIAL  DISTRIBUTION  LIST 

Addressee  No.  of  Copies 

ASN  ( RE&S)  I 

OUSDR&E  (Research  and  Advanced  Technology)  1 

Deputy  U5DR&E  (Res  &  Adv  Tech)  1 

OASN ,  Spec  Dep  for  Adv  Concept  1 

Principal  Dep  Assist  Secretary  (Research)  1 

OASN,  Dep  Assist  Secretary  (Res  &  Applied  Space  Tech)  1 

ONR,  0NR-100,  -102,  -200,  -400,  -410,  -411  (N.  Gerr) ,  9 

-422,  -425AC ,  -430 

CNO,  OP-098 ,  -941 ,  -951  3 

CNM,  MAT-00,  -05,  -05  (G.  Morton)  3 

DIA,  DT-2C  1 

NRL,  Code  5162  (Nate  Yen),  (A.  A.  Gerlach)  2 

NRL,  USRD  1 

NORDA  1 

USOC,  Code  241,  240  2 

SUBBASE  LANT  1 

NAVSUBSUPACNLON  1 

NAVOCEANO,  Code  02  (1)  2 

NAVAIRSYSCOM  1 

NAVELECSYSCOM,  ELEX  03,  310,  320  3 

NAVSEASYSCOM,  SEA-00,  -05R,  -06R,  -63R,  SEA-92R  5 

NAVAIRDEVCEN,  Warminster  1 

NAVAIRDEVCEN ,  Key  West  1 

NOSC,  Code  713  (Fredric  J.  Harris),  8302  2 

NOSC ,  Code  6565  (Library)  1 

NAVWPNSCEN  1 

NCSC ,  Code  724  1 

NAVCIVENGRLAB  1 

NAVSWC  1 

NAVSURFWPNCEN ,  Code  U31  1 

NISC  1 

CNET,  Code  017  1 

CNTT  1 

NAVSUBSCOL  1 

NAVTRAEQUIPCENT,  Technical  Library  1 

NAVPGSCOL  1 

NAVWARCOL  1 

NETC  1 

APL/UW,  SEATTLE  1 

ARL/PENN  STATE,  STATE  COLLEGE  (F.  W.  Symons)  1 

CENTER  FOR  NAVAL  ANALYSES  (ACQUISITION  UNIT)  1 

DTIC  1 

DARPA  1 

NOAA/ERL  1 

NATIONAL  RESEARCH  COUNCIL  1 

WEAPON  SYSTEM  EVALUATION  GROUP  1 

WOODS  HOLE  OCEANOGRAPHIC  INSTITUTION  (Dr.  E.  Weinstein)  1 

ENGINEERING  SOCIETIES  LIB,  UNITED  ENGRG  CTR  1 

NATIONAL  INSTITUTE  OF  HEALTH  1 


TR  7377 


INITIAL  DISTRIBUTION  LIST  (Cont'd) 

Addressee  No.  of  Copies 

ARL,  UNIV  OF  TEXAS  1 

MARINE  PHYSICAL  LAB,  SCRIPPS  1 

UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO  1 

NAVSURWEACTR  1 

DELSI  1 

INTERMS  INC.  1 

MAR  INC.  1 

B-K  DYN  INC.  1 

BBN,  CAMBRIDGE,  MA  (H.  Gish)  1 

EWASCTRI  1 

HYDROINC  (D.  Clark)  1 

SUMRESCR  (M.  Henry)  1 

ANALTECH  INC.  (G.  Bourgond)  1 

ANALTECHNS  1 

GENPHYCORP  (M.  Bauer)  1 

EDOCORP  (J.  Vincenzo)  1 

OPERRES  INC.  (Dr.  V,  P.  Simmons)  1 

TRA  CORP.  (J.  Wilkinson)  1 

BB&N  INC.,  San  Diego  1 

NETS  (R.  Medeiros)  1 

GESY  (D.  Bates)  1 

ORI  CO.,  INC.  (G.  Assard)  1 

Dr.  Norman  Bleistein,  Denver,  CO  1 

PROMETHEUS  INC.  (Dr.  J.  S.  Byrnes)  1 

BB&N  INC.,  New  London,  CT  (Dr  P  G.  Cable)  1 

ROYAL  MILITARY  COLLEGE,  Dept,  of  Electrical  Engineering, 

Kingston,  Ontario,  Canada  (Prof  Y.  T  Chan)  1 

UNIVERSITY  OF  FLORIDA,  Dept,  of  Electrical  Engineering 

(D.  C.  Childers)  1 

SANDIA  NATIONAL  LAB,  Albuquerque,  NM  (J.  Claasen)  1 

COGENT  SYSTEMS,  INC.,  Waban,  MA  (J.  P.  Costas)  1 

BB&N  INC.,  Arlington,  VA  (Dr.  H.  Cox)  1 

AUTONETICS  MARINE  SYSTEMS  DIV.,  Anaheim,  CA  (L.  T.  Einstein)  1 

ROCKWELL  INTERNATIONAL,  Anaheim,  CA  (Dr.  D.  F  Elliot)  1 

Bernard  Harris,  Dobbs  Ferry,  NY  1 

EG&G,  Manassas,  VA  (Dr.  J.  Hughen)  1 

MAGNAVOX  GOV.  &  IND.  ELEC.  CO.,  Dept.  525,  Fort  Wayne,  IN 

(R.  Kenefic)  1 

DREXEL  UNIVERSITY,  DEPT.  OF  ELEC.  &  COMP.  ENG., 

Philadelphia,  PA  (Prof.  S.  Kesler)  1 

DEFENCE  RESEARCH  ESTABL.  ATL.,  Dartmouth,  Nova  Scotia, 

Canada,  Library  (B.  E.  Mackey)  1 

SCHLUMBERGER  WELL  SERVICES  ENGINEERING,  Houston,  TX 

(S.  L.  Marple)  1 

PSI  MARINE  SCIENCES,  New  London,  CT  (Dr.  R.  Mellen)  1 

Dr.  David  Middleton,  New  York,  NY  1 

UNIVERSITY  OF  CONNECTICUT,  Dept,  of  Elec.  &  Comp.  Science 
(Prof.  C.  L.  Nikias) 


1 


TR  7377 


INITIAL  DISTRIBUTION  LIST  (Cont'd) 


Addressee  No.  of  Copies 

BB&N  INC.,  Canoga  Park,  CA  (Dr.  A.  G.  Piersol)  1 

LINK-A-BIT  CORP.,  Lexington,  MA  (Dr.  R.  Price)  1 

DALHOUSIE  UNIVERSITY,  Dept,  of  Oceanography,  Halifax, 

Nova  Scotia,  Canada  (Dr.  B.  Ruddick)  1 

UNIVERSITY  OF  RHODE  ISLAND,  Dept,  of  Electrical  Engineering 

(Prof.  L.  Scharf,  Prof.  D.  Tufts)  2 

YALE  UNIVERSITY,  Depth  of  Electrical  Engineering 

(Prof.  P.  M.  Schultheiss)  1 

DEFENSE  SYSTEMS,  INC.,  McLean,  VA  (Dr.  G.  S.  Sebestyen)  1 

NATO  SACLANT  ASW  RESEARCH  CENTRE  (Dr.  E.  J.  Sullivan,  Library)  2 

DEFENCE  RESEARCH  EST.  PACIFIC,  Canada  (Dr.  D.  J .  Thomson)  1 

Robert  J.  Urick,  Silver  Spring,  MD  1 

RCA/GOVERNMENT  SYSTEMS  DIV.,  Moorestown,  NJ  (H.  Urkowitz)  1 

TEL-AVIV  UNIVERSITY,  Dept,  of  Electronic  Systems,  Tel-Aviv, 

Israel  (Prof.  E.  Weinstein)  1 

TRACOR,  Inc.,  Analysis  &  Appl .  Res.  Div.,  Austin,  TX 

(J.  F.  Wilkinson)  1 

COAST  GUARD  ACADEMY,  Dept,  of  Mathematics  (Dr.  J.  J.  Wolcin)  1 

UNIVERSITY  OF  IOWA,  Elec.  &  Comp.  Eng.  Dept. 

(Prof.  D.  Youn)  1 


U219560 


I 


4 


